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MONTE CARLO METHOD FOR THE CALCULATION OF TRANSPORT 
PROPERTIES IN A LOW-DENSITY IONIZED GAS 
by Charles M. Goldstein 
Lewis Research Center 

SUMMARY 


An introduction to the general Monte Carlo method is presented, along with a dis- 
cussion of its scope of application to plasma physics. This is followed with a heuristic 
sketch of the method. The problem of electron flow through a perfect Lorentzian gas in 
a parallel- plane diode is then formulated. The Monte Carlo solution is discussed in de- 
tail along with the relevant computational techniques employed. (Pertinent concepts 
from the theory of random variables are included in an appendix. ) 

The effects of mean free path on current-voltage characteristics, density distribu- 
tion, and potential distribution are presented for two cases - monoenergetic and thermi- 
onic emission. Results indicate that electron-neutral elastic collisions can have a sig- 
nificant effect on the current-voltage characteristics for electrode separations as small 
as one mean free path in the case of thermionic emission, and one-half mean free path 
in the case of monoenergetic emission. 


INTRODUCTION 

A major difficulty in the study of low-density ionized gases is the lack of suitable 
analytical methods for determining the effects of collisions. "Low density" is here 
defined as those situations in which a characteristic dimension is of the order of a few 
mean free paths, that is, the regime wherein neither a collisionless nor a continuium 
approximation can be expected to represent the actual situation. This regime is of im- 
portance in low-pressure thermionic diodes, plasma sheaths (probe theory), ion engines, 
and cross-section measurements. 

Much effort has been expended to obtain solutions of the Boltzmann transport 
equation for low-density neutral gases (ref. 1). Little work has been done on the 
extension of these methods to low-density ionized gases. Recently, however, Sockol 



(ref. 2) has succeeded in numerically integrating the Boltzmann transport equation for 
a particular low-density ionized gas problem. Unfortunately, the numerical integration 
is very difficult for even the simplest hard-sphere collision cross section; the feasi- 
bility of extending this method for more complex collision cross sections has not as yet 
been investigated. 

This report presents a new method for determining analytically the transport prop- 
erties in a low-density ionized gas for an arbitrary collision cross section. Results 
with this method are given for two electron transport problems. The method proposed 
is, essentially, a consistent-field Monte Carlo method. Since the Monte Carlo method 
has not been widely applied in the fields of plasma physics or ionized gases, a brief 
introduction is presented; a review of the pertinent random variable theory is given 
in appendix A. The ability to use this method effectively is strongly dependent on 
numerical procedures and ’’tricks of the trade.” A section is therefore included on 
the various computing techniques; in addition, a complete program listing plus selected 
flow charts are to be found in appendix C. A discussion of two important computing 
programs is presented by their author, H. Renkel, in appendix B. 

The general method of calculating transport properties is applied herein to the 
problem of electron transport in a perfect Lorentzian gas with a hard- sphere collision 
cross section; in particular, the method is employed to obtain electron flux character- 
istics in a plane -parallel diode including the effect of electron-neutral collisions. 

Langmuir (ref. 3) published the first correct solution to the effect of space charge 
and initial velocities on the potential distribution and thermionic current between 
parallel -plane electrodes for no collisions (vacuum diode). He also studied (ref. 4) 
the problem of diffusion of electrons back to the emitter for the case of a very small 
mean free path. These results are extended herein to the case of electron- neutral 
collisions for which the mean free path is not necessarily small with respect to the 
inter electrode separation. 


MONTE CARLO METHOD 
General Description 

The Monte Carlo method is, in general terms, a technique for solving physical and 
mathematical problems by using random sampling. Although the term ’’Monte Carlo 
method” has been subjected to various interpretations, an acceptable statement of the 
method as applied herein has been given by Donsker and Kac (ref. 5): ’’The Monte Carlo 
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approach consists in permitting a ’particle* to play a game of chance, the rules of the 
game being such that the actual deterministic and random features of the physical proc- 
ess are step by step exactly imitated by the game. By considering very large numbers 
of particles, one can answer such questions as the distribution of the particles at the end 
of a certain period of time, the number of particles to escape through a shield of a spe- 
cific thickness, etc. One important characteristic of the preceding approach is that the 
functional equation describing the diffusion process is bypassed completely, the proba- 
bility model used being derived from the process itself. ” 

A short history of Monte Carlo applications is to be found, in the paper by Goertzel 
and Kalos (ref. 6). An excellent review of the basic principles is given in reference 7, 
and an extensive bibliography has recently been compiled by Kraft and Wensrich (ref. 8). 
This method has, in recent years, been employed with considerable success to a wide 
variety of problems, most notably in the area of nuclear shielding problems (viz. , neu- 
tron transport). These latter problems are linear in the sense that the neutron trajec- 
tories are independent of the neutron density. More recently, the method has been ex- 
tended to certain nonlinear problems in radiation transport (ref. 9). 

There have been some applications of the Monte Carlo methods to the investigation of 
one -dimensional electron (ion) diodes, but these studies are more often referred to as 
computer-simulated solutions, or ’’computer experiments.” The difference in termi- 
nology reflects the fact that these studies approximate the physical model by a finite 
number of current sheets, which are then followed deterministically through all mutual 
interactions by the computer. The Monte Carlo method, on the other hand, most fre- 
quently implies repeated, stochastically independent trials. A short history of the afore- 
mentioned computer experiments is to be found in the paper by Burger (ref. 10). These 
studies do not take collisions into account, nor does it seem practical to do so because 
of the demands this type of analysis would impose on computer storage requirements. 

Itoh and Musha (ref. 11) employed a Monte Carlo calculation to determine the ion- 
ization and excitation coefficients of electrons in a uniform electric field E for given 
gas pressure P. They also computed drift velocity and mean energy for several values 
of E/P. Although the authors state that this method can be extended to strong, non- 
uniform electric fields, it cannot provide a suitable model from which diode character- 
istics could be obtained since space-charge effects, which introduce a nonlinearity, have 
not been considered. 

Just as the nonlinearity in the radiation transport problem is characterized by a 
single parameter, the temperature (ref. 9), so the nonlinearity in the charge-particle 
transport problems is characterized by the potential. Unlike the photons in the former 
problems, however, the charged particles experience a body force proportional to the 
first derivative of the potential. 


k. 
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Heuristic Sketch of the Method 


First the relatively simple problem of the attenuation of a molecular beam by a 
homogeneous gas shall be considered. If the actual experiment is performed for a given 
emission flux r Q and the flux reaching the target (or collector) r c is measured, it is 
reasonable to interpret the ratio r„/r_ as the probability that a unit of emitted flux 
will reach the target. If a knowledge of the scattering probabilities of a single molecule 
passing through the same gas is assumed, it is possible on a computer to follow a cer- 
tain number N , one at a time, and tally the number N„ that reach the collector (the 
others are scattered back to the emitter). Then the ratio N /N would be an approxi- 
mation to the experimentally determined r c /r Q . Since the experimental fluxes may be 
of the order of 10 ^ particles per square centimeter per second or higher, it is not con- 
ceivable, even in this relatively simple situation, to "do the experiment” on the com- 
puter. As N q becomes larger, however, the approximation N c /N q becomes better. 
Statistical analysis provides a means of estimating how good the approximation is. 

For the case of charged particles flowing through a gas, the situation is complicated 
by the nonlinearity introduced by the space charge. That is, the flow of charged par- 
ticles is not only influenced by collisions with the gas molecules, but also by the poten- 
tial field; the potential field, itself, is a function of the density of charged particles. 

This is the situation considered herein. To start, for a given collector potential T(L) 
a potential distribution is assumed. An approximation to the current reaching the 

collector, N./N., is then obtained as in the molecular beam case. In addition, however, 
the contribution to the density, at preselected data points, of each charged particle is 
also tallied. These densities are used to solve Poisson’s equation for a new potential 
distribution. The process is then repeated (i. e. , iteration is performed on the potential 
distribution) until the potential distribution "converges. ” Convergence must here be 
considered only in a statistical sense; when further iterations produce only random 
fluctuations in the potential distribution, "convergence” is assumed. Random fluctua- 
tions are, of course, to be expected, since only a small statistical sample N Q of the 
total flux is considered. After convergence has been achieved, succeeding iterations 
may be considered, in the parlance of statistical analysis, as independent trials; each 
resulting approximation N /N may be considered a sample mean (appendix A). An 
analysis of the sample means provides a way of estimating the accuracy of the result. 


FORMULATION OF ELECTRON DIODE PROBLEM 

The physical model of an infinite parallel -plate diode is depicted in figure 1. In the 
same figure the types of scattering that may occur for a monotonic potential distribution 
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are shown. In figure 2 a typical potential 
distribution is shown. When a potential 
minimum exists as indicated in this figure, 
a certain portion of thermionically emitted 
electrons will be rejected back to the 
emitter even in the absence of collisions. 
The existence of a potential minimum less 
than both emitter and collector potentials 
defines the space-charge-limited regime 
of diode operation. 

The perfect Lorentzian gas assump- 
tion implies an infinite mass target par- 
ticle, and hence the laboratory system be- 
spatiai coordinate, x comes equivalent to the center of mass 

Figure 1. - Diode model and types of scatter. system. Since hard- sphere collisions 

result in isotropic scattering in the center 
Collector Q f m ass system, the equivalence of the 
■i (L) two systems in this case results in iso- 
tropic scattering in the laboratory system. 

Isotropic scattering means, by defini- 
tion, that the probability of scattering into 
unit solid angle is the same for all angles. 
The probability distribution function 

0 x minimum L (hereafter, p. d. f. , see appendix A) is 

Spatial coordinate, x 

therefore the constant 1/4jt. Hence the 

Figure 2. - General potential distribution. 

probability of scattering into solid angle 
dS2 is dS2/47T. In terms of the scattering angle 0, this becomes 

™ = (i) 

2 

Consequently, the p. d. f. of scattering into angle d0 is simply sin d/2. 

The assumption of hard-sphere collisions also implies a constant mean free path X. 
Now if that group of electrons that has just collided is considered, then the fraction of 
these electrons that will suffer collisions in distance i is (ref. 12, p. 102, eq. (98a)) 

1 - e~ aS -c (2) 

where 
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a = s/l 


a = L/x 


( 3 ) 


and s is the path length and L the electrode spacing. 


Thermionic Emission 

For this case it is assumed that electrons are emitted with a half-Maxwellian 
velocity distribution 


2 2 

f(u, V) du dV = (4/^7T)Ve"^ u + V ^ 


du dV 


where 


u = v x / y2kT/ m 




V = 


y z 

yikrTm 


The one -dimensional Poisson’s equation, in dimensionless variables, becomes 


(4) 


(5) 


<P"{ y) = C n(y) 


where 


and 



_ eV 
<P= — 

kT 


> 



n o J 


( 6 ) 


( 7 ) 


6 



C = 8(7r/2kT) 3y/2 m 1//2 eJ 0 L 2 


( 8 ) 


Monoenergetic Emission 

The analysis for monoenergetic emission directly parallels that for thermionic col- 
lision with a few minor changes. The dimensionless variables u, V, and <p are now 
defined as 


u = 


v.VS 


<P = 


2eV 


mv. 


(9) 


The constant parameter C in Poisson’s equation (eq. (6)) becomes 


8jreL 2 J 


( 10 ) 


mv. 


MONTE CARLO SOLUTION 
Thermionic Emission 

Initial conditions. - It must be emphasized that the test ’’electrons” are not chosen 
from the half- Maxwellian distribution (eq. (4)). Although test ’’electrons” are men- 
tioned, the statistics are obtained for units of electron flux - not units of charge. Hence, 
the initial velocities must be chosen from the distribution of flux in velocity space 

f* uf(u, V) du dV = 4uVe-(“ 2+v2 ) d “ dv <“> 
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from equation (4). Since the u and V components of velocity are independent, the 
respective marginal distributions (see ref. 13, p. 287) F u (V) and Fy(u) can be 
obtained: 


T (V) m f du f V 4uVe"^ u2+V ^ dV 


( 12 ) 


r°° r u 
r(u) = / d V 

J 0 J 0 


4uVe 


-(u 2 +V 2 ) 


du 


(13) 


But these marginal distributions are simply the cumulative distribution functions (here- 
inafter, c.d. f. ) for u and V, respectively. From equations (A 16) and (A15), 


u 2 = -ln(R u ) 
V 2 = -ln(Ry) 


(14) 


where R u and Ry are random numbers between 0 and 1. Equations (14) are then used 
to determine the initial velocities of each test electron. 

Distance to collisions . - The distance to collision must be obtained at the start of 
each new electron trajectory (i. e. , on emission from emitter or after a collision). 

Equation (2) can also be interpreted as the probability that an electron will suffer a 
collision in a distance f < j£ . This, however, is just the definition of the c.d. f. F(£ ) 
(see appendix A). Hence, from equation (A16) can be obtained a relation between the 
random numbers R^ and the distribution of path length to collision: 

-al 

R f = 1 - e c (15) 


where 


tc = .(i) ^ . Rl) 
= - 0 ln( ^ 


(16) 


Scattering angle. - If a collision takes place, then the scattering angle 9 


must be 
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determined. From the p. d. f. (eq. (1)), thec.d.f. F(0) can immediately be obtained 
(compare eq. (A3) in appendix A): 


F(0) = — - cos - ? 
2 


( 17 ) 


But in this case the c. d. f. can take on the values -1 < F(0) < 1 (forward and backward 
scattering). Hence, in order to choose randomly from this range (see eq. (A16)), let 


cos 0 = 1- 2R 0 


(18) 


2 

where once again 0 < < 1. Equation (18) is the final result since only cos 0 (and 

2 2 y 

sin 0 = 1 - cos 0) is of interest in the actual computations. 

Charge density . - The data points y i are selected by the curve -fitting subroutine 
(see appendix B where y^ corresponds to the arguments x a ). The contribution of the 
k th test electron (unit of flux) of velocity u k (y) to the charge density at each is 


n k (Yi) = 



(19) 


where 


u(y) = 



( 20 ) 


y Q is the position of the last "event" (collision or emission), and u Q is the initial 
velocity immediately after the last event (i. e. , at the beginning of a new trajectory). 
The tallied density at a data point y^^ for a total of N q histories is then 


n( Yi ) = 



1 

u k^i> 


( 21 ) 


where the sum over k may be greater than, equal to, or less than N Q because of col- 
lisions and turning points in the potential field. 

Current to col lector . - The ratio of current density to the collector J to the 
emitted current density J Q for each iteration is computed from the relation 
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( 22 ) 



where N c is the number of test electrons reaching the collector. 


Monoenergetic Emission 


The solution for monoenergetic emission is exactly the same as for thermionic 
emission with two exceptions. First, instead of choosing from an initial distribution of 
velocities the initial conditions are, for every electron, u = 1 and V = o, and second, 
the density for N Q histories becomes (cf. eq. (14)) 



(23) 


COMPUTATIONAL METHODS 

The original program was based on the assumption that Monte Carlo computations 
would be limited to only a few collisions because of the requirement of reasonable com- 
putor execution times. Hence, this program was optimized for L/A < 1. The results 
proved this assumption overly pessimistic, but the program was not revised for the 
present report. 


Evaluation of Potential and Potential Minimum 

After a curve fit of the density is obtained (subroutine CHEBY, appendix B), the 
density distribution is given by a power series in y: 

n(y) = a Q + a x y + a 2 y 2 + . . . + a n y n 

k 

= Y, a / J (24) 

j=0 

where k is the degree of the fit. 
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The potential is obtained by substituting equation (24) in Poisson’s equation (eq. (6)) 
and integrating n(y) term by term: 


cp(y) = c Q + cy + 


CV h yj-2 

j + D(j + 2) 


(25) 


But since 


<p(o) = 0 


equation (25) can be written 


k+2 

<p{ y) =^ c jY j ( 26 ) 

j=l 



a j - 2 
j(j - 1) 


i > 2 


k+2 

Cj = <p(l) -^^ C J 
j=2 


-S 




(27) 


After the a^ have been determined (subroutine CHEBY), the Cj are computed in sub- 
routine COEF(2). 

Originally, equation (26) was employed (with k usually equal to 10) each time cp( y) 
was evaluated, but this proved too time consuming. For this reason it was decided to 
tabulate <p(y) at the beginning of each iteration and use the tabulated values whenever 
possible. The interelectrode space was subdivided into 1024 regions, and the 1025 val- 
ues of cp( y) were tabulated in subroutine MINPHI. At the same time, cp was tested at 
each evaluation for the minimum value. Hence, the location of the potential minimum 
was ascertained within ±1/2048 of the interelectrode separation. In addition, cp was 
tabulated at data points y^ where the density was to be tallied. The results of tabu- 
lating the potential was an eight-fold (and greater) decrease of execution time. 
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Choosing from the Distribution e~ x dx 


It is pointed out in appendix A that choosing random values X k from the distribu- 
tion whose p. d. f. is e -x dx is equivalent to choosing random numbers R k from the 
uniform distribution (eq. (A7)) and using equation (A15): 

X k = -ln(R k ) (28) 

2 2 

In the present problem, it is possible to identify the random variables u and V 
with X, and £ with (1/ «)X (eqs. (14) and (16), respectively). The random num- 

L 

bers R k are obtained from a pseudo-random-number generator of the congruence- 
method type (ref. 5). This random generator is part of the computor library here at 
Lewis Research Center. 

Although the desired random variable can be obtained directly from equation (28), 
it was decided to tabulate the X k instead. A table (1025 entries) was constructed of 
X k (subroutine CUMVEL) at the beginning of the program. The table look-up is five 
times as fast as employing equation (28) each time. 


Location of Collision 

If a distance to collision is given, the location of the collision y c is obtained 
by solving 


^c 


r 

•v \ % + rfy 


dy' 


(29) 


) - <p(y Q ) 


Two methods were used to minimize the number of times the integrand, and specifically 
cp{ y), need be evaluated. First, Simpson's rule was used in a search routine to allow 
the use of the tabulated values of cp( y) (see Evaluation of Potential and Potential Mini- 
mum section, p. 10), and then the step size (in the use of Simpson's rule) was made to 

Op ' 

depend on the ratio V /u . 

The procedure employed for obtaining a reasonable step size can be best explained 
by an example. Assume that y c falls between any two points y Q and y f . For a 
straightforward application of Simpson's rule, three values are needed of the integrand 
in equation (29) at three equidistant values of y: y 1? y 2 , y g . Initially y^^ = y The 
<p( y) has already been tabulated at 1025 values of y given by 
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m 

1024 


m = 0, 1, 2 


1024 


(30) 


Consequently, y Q and y^ will always be selected equal to tabulated values of y m and 
y . Hence, the first estimate of step size A in units of m is given by 0 


A = 


m f - m 
f o 


(31) 


where [ ] refers to the integral value. A second estimate of step size (obtained as a 
result of trial and error computations) is given by 


A’ = 2 


4-M 


where 


M = [log 10 (V 2 /u 2 )] 


Then, the step size is taken as the minimum of the two estimates. 

If the value of A from equation (31) is zero (i. e. , distance to collision is less than 
four steps), then A is set equal to 


A = 


m f - m o 


If this should be zero, then the collision location y is arbitrarily set equal to y 


m x 


CONVERGENCE AND STANDARD DEVIATION 

It was observed, during tests of the program, that convergence (in the statistical 
sense, p. 4) was obtained in the first few iterations. Since the succeeding iterations are 
treated as independent trials, the problem arises in a production run of just how to 
decide when convergence occurs. This was done in the following manner. 

Each case (given anode potential) was run for a given number of iterations, for 
example, 15. At the end of each iteration the sample means n(y.) and J/J Q (see 
eqs. (21) to (23)) were stored. Each of these stored values is analogous to an experi- 
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TABLE I. - EFFECT OF VARIOUS PARAMETERS ON STANDARD DEVIATION AND EXECUTION TIME 


Item 

Type of 
emission 

Electrode 
spacing 
to mean 
free path 
ratio, 
L/A 

Dimension- 
less collec- 
tor potential, 

*(i> 

Current 

density 

ratio, 

j/j 0 

Standard 

deviation, 

a j 

Sample 

size 

Number of 
iterations 

Collisions 
(for one 
iteration) 

Execution 

time, 

min 

1 

Mono- 

0. 

1 

0. 75 

0. 961 

0.001 

5 000 

10 

483 

2.42 

2 

ener- 



2.0 

. 972 

.0012 

5 000 

10 

503 

2.46 

3 

ge 

tic 



2.0 

. 971 

.0019 

2 000 

5 

213 

.51 

4 





4.0 

. 986 

.0005 

10 000 

10 

1 079 

5.28 

5 




5 

4.0 

. 918 

.0015 

2 000 

10 

1 257 

2.48 

6 



1 

0 

4.0 

. 839 

.0015 

1 000 

10 

1 483 

2.09 

7 

The 

"mi- 

5 

0 

12.0 

. 358 

.009 

1 000 

15 

8 475 

17. 7 

8 

onic 

5 

0 

32.0 

.672 

.003 

1 000 

18 

16 367 

37.65 

9 




1 

10.2 

. 942 

.0016 

1 000 

10 

147 

1.49 

10 




1 

10.2 

. 942 

.0008 

10 000 

10 

1 198 

14. 18 


mental data point. Carrying the analogy further, at the end of 15 runs (iterations) there 
were sets of 15 data points for each of the sample means. If this were an experiment, 
it would be expected that each set of 15 data points would have a certain amount of 
’’scatter” due to random error. In the present situation, however, the iterations before 
convergence will produce data points with a nonrandom error. The problem then be- 
comes one of simply eliminating the iterations that introduce a nonrandom error. This 
was accomplished by obtaining the sample mean and standard deviation (see appendix A) 
of each set of 15 data points. Then from each set only those points were retained that 
were within three standard deviations of the sample mean. The final values of sample 
means and standard deviation (given in table I) were obtained from the remaining data 
points. In all cases, the number of iterations treated as independent trials was of the 
order of ten. 


RESULTS 

Thermionic Emission 

The effect of mean free path on the current-voltage characteristic is shown in fig- 
ure 3. The solid line, L/A = 0, represents the collisionless solution of Langmuir 
(ref. 3). The Monte Carlo calculations indicated along this curve were undertaken as a 
check on the computer program. These particular results were obtained with 5000 his- 
tories per iteration and ten iterations. The execution time for each point on the curve 
varied between 2. 5 and 4. 0 minutes. 
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Figure 3. - Effect of mean free path on current-voltage characteristics for thermionic emission. Dimensionless 
constant C * 50. 


The two solid data points on the curves for L/A = 1 and 5 represent the conditions 
where the slope of the potential is zero at the emitter. The O’s on the curve L/A = 1 
indicate the results of an independent solution of Boltzmann’s transport equation for this 
problem (ref. 2). 

The curve for L/A = 5 was not extended to lower cp( 1) because of a loss in pre- 
cision in the curve-fitting routine program (appendix B) used to fit the density distribu- 
tion. A more flexible routine is being developed. 

The effect of potential on the electron density distribution is shown in figure 4. 

From the emitter out to about one mean free path, the density of the higher energy elec- 
trons is less than that of the lower energy electrons as would be expected under condi- 
tions of no collisions. The actual decrease in the magnitude of the density at the emitter 
surface, however, indicates that in the higher potential case more of the backscattered 
electrons are being turned about by the potential field before reaching the emitter. This 
can be best understood by considering the effect of an accelerating potential field on the 
cone of capture at the emitter for backward scattering (see fig. 5). This cone of capture 

may be defined by a polar angle 6*. It will suffice to consider a first collision whereby 

2 2 

the electron has initial energy of ui| + and the collision occurs at x . The magni- 

O O Cr 

tude of the x-component of velocity after scatter becomes 

u 2 = [ u o + V o + ^ x c>] cos2(9 ( 32 ) 
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Figure 4. - Effect of anode potential on electron density distribution for therm- 
ionic emission. Dimensionless constant C = 50; electrode spacing to mean 
free path ratio LM = 5. 


If the electron is to reach the emitter against the monotonic potential field <p(x) > 0, 
2 

then u must satisfy the condition 


u 


2 


^ <?(x c ) 


(33) 


When equation (33) is substituted into equation (32), 9* is defined by 


cos 9* = 


^ x c ) 

u 0 + V o + ^ x c) 


Emitter Collector 



Figure 5. - Cone of capture at emitter. 


or 


cos 2 0* = (34) 

u 2 + V 2 
1 + — ° 

cp(x c ) 

Equation (34) shows that an increase in poten- 

o 

tial <p(x c ) increases cos 9* and reduces 9*. 
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Figure 6. - Effect of mean free path on electron density distribution for thermionic emission. Dimensionless con- 
constant C = 50; dimensionless collector potential <p( 1) = 32. 


Thus, the higher the potential field, the smaller is the cone of capture at the emitter. 

This same phenomenon accounts for the crossover in the curves of figure 4 away 
from the emitter. Since fewer of the backscattered electrons in the higher electric field 
case reach the emitter, this implies that more are turned about by the field. The pres- 
ence of turning points in the electron trajectories affects the charge density in two ways. 
Since the u-component of velocity is zero at a turning point, the contribution to the 
charge density there is exceptionally high; in addition, the path length of an electron in 
the neighborhood of a turning point is much greater than the distance traveled normal to 
an electrode surface, these electrons suffer more collisions, and, hence, contribute 
more strongly to the charge density. This latter point is vividly illustrated by compar- 
ing the typical number of collisions per iteration for the two cases of figure 4 (items 7 
and 8, table I, p. 14). In the low potential case (<p( 1) = 12) over 8000 collisions were 
observed in one iteration, while for the high potential case (<p(l) = 32) over 16 000 col- 
lisions were observed. The increase in number of collisions accounts for the crossover 
in the two curves of figure 4 and the higher density for y > 0. 2 in the case <p{ 1) = 32. 

The effect of mean free path on the density and potential distributions for constant 
collector potential are shown in figures 6 and 7, respectively. As expected, the effect 
of collisions is to increase the charge density and, therefore, decrease the potential in 
the interelectrode space. 
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Figure 7. - Effect of mean free path on potential distribu- 
tion for thermionic emission. Dimensionless constant 
C = 50. 


Monoenergetic Emission 

The corresponding diode characteristics 
for monoenergetic emission are shown in fig- 
ures 8 to 11. The author has, at present, no 
hypothesis regarding the inflections observed 
in the current-voltage characteristics (fig. 8) 
for Li/A = 0. 5 and 1.0. The points calculated 
are reproducible, and each point, as plotted, 
spans at least plus or minus two standard 
deviations about the mean J/J Q . The solid 
lines represent independent solutions of the 
Boltzmann equation for this problem in the 
limit of one collision (ref. 14). 

Another noteworthy feature of the mono- 
energetic emission characteristics is the 
buildup of charge density in the interelectrode 
region as the potential is decreased (fig. 9). 
This increase in charge density is consider- 
ably enhanced by the appearance of a potential 
minimum (upper curve in fig. 9). The poten- 
tial minimum causes more turning points to 
occur in the trajectories of the scattered 



Figure 8. - Effect of mean free path on current-voltage characteristics for monoenergetic 
emission. Dimensionless constant C = 10. 
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Dimensionless electron density distribution, niy) Dimensionless electron density distribution, n(y) 





electrons. Since the u -component of velocity 
becomes zero at a turning point, the contribu- 
tion to the charge density of electrons under- 
going reflections in the potential field is ex- 
ceptionally high. 

DISCUSSION OF RESULTS 

The agreement of the solution obtained by 
the method proposed in this report and the in- 
dependent results of Sockol (ref. 2 and fig. 3, 
p. 15) and Goldstein and Goldstein (ref. 14 and 
fig. 8) is very gratifying indeed. Most en- 
couraging, with regard to the extension of this 
method to other problems, are the statistics 
presented in table I (p. 14). These statistics 
show that the execution times needed to obtain reasonable standard deviations Oj need 
not be excessive. In turn, they illustrate the effect of the consistent -field constraint 
(Poisson’s equation) on the number of histories needed for good statistics (tens and hun- 
dreds of thousands of histories are generally required in other problems where this con- 
straint is absent). It must be emphasized that the execution times illustrated in table I 
are not the minimum attainable, since no attempt has yet been made to incorporate any 
of the variance -reducing techniques discussed in the literature (ref. 15). 

The execution times in the problems treated herein could most directly be de- 
creased by a more extensive use of tabulated values (eliminating the Gaussian quadra- 
tures - hence, obviating completely the need to evaluate </?(y) during an iteration) and 
by optimizing the number of tabulations needed (one may not need 1025 tabular values). 

In addition, for larger values of L/A, it would be more appropriate to step along each 
trajectory from the emitter instead of first ascertaining if a collision has occurred in 
the interelectrode space as is done in the present case. 


Figure 11. - Effect of mean free path on potential distribu- 
ion for monoenergetic emission. Dimensionless constant 
C ■ 10/ Vi. 


CONCLUDING REMARKS 

A general method for the calculation of transport properties in a low-density ionized 
gas has been presented. This method has been applied to two cases of electron trans - 
port in a perfect Lorentzian gas. Excellent agreement has been demonstrated by two 
other independent investigations. 

Although the particular applications of the method presented herein employ a hard- 
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sphere collision model, the great advantage in the Monte Carlo method lies in its in- 
herent ability to provide similar solutions for any given collision model, theoretical or 
experimental. This includes inelastic, charge exchange, and ionizing collisions. This 
method is limited, however, to those cases where avalanche ionization does not occur; 
even in this latter case, however, the Monte Carlo method should be capable of provid- 
ing the source intensities for the collision term in the Boltzmann equation for arbitrary 
cross sections, and, therefore, allow a numerical solution of the same. 

This method should also be of value in the solution of plasma sheath problems, 
which are, in reality, just generalizations of the diode problem with different boundary 
conditions at the emitter and/or collector. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, May 14, 1965. 
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APPENDIX A 


RANDOM VARIABLES 

In this section there is interest only in continuous probabilities for which there 
exists a continuous function f(x) , called the probability distribution function (hereinafter 
p.d. f. ), such that 


P[a < X < 



/ +°o 

f(x)dx = 1 

on 


(Al) 


where 


f(x) >0 -°° < X < oo 

To speak of a "random variable" X (instead of x) is really to define a mathe- 
matical point of view. This unambiguous point of view maintains no interest in the exact 
value of X but instead is only interested in inquiring about the probability of finding X 
in a certain region (of x-space). 

For example, the case is considered where the probability density function is the 
nondimensionalized Maxwellian distribution of the x-component of flux (eq. (11) inte- 
grated over V): 

-u 2 

f(u)du = 2ue du u > 0 

(A2) 

= 0 u < 0 

In the present analysis, the concern is not for a knowledge of a particular value of u, 
but rather to determine just what the probability is that a random variable U lies in the 
range u, u + Au. 

In the study of a random variable X the function F^(x) is of great importance: 



This F-j^(x) , or simply F(x), is called the cumulative distribution function (hereinafter 
c. d. f. ) of the random variable X. This function shall be used subsequently. 

A concept basic to the discussion of random variables is the expectation value E[ ] 
of a function of a random variable g(X): 

/ +°° 

g(x)f(x)dx (A4) 

OO 

In particular, the expectation of a random variable itself 

/ +°o 

xf(x)dx (A5) 

■OO 


is the familiar mean or average value of X, -«> < X < Of interest in the text is the 
expectation of the function l/U of the random variable U distributed as f(u) (eq. (A2)). 


hence, 



(A6) 


Choosing from a Distribution 

It is first necessary to define what is meant by choosing a sequence of random num- 
bers X^ from a distribution f(x) (or equivalently, choosing X^ distributed as f(x)). 

It is assumed, for the sake of illustration, that the p. d. f. f(x) is nonzero only in the 
interval Osxsl. This interval is then subdivided into 10 equal subintervals. Then, 
if the sequence of N random numbers X t is distributed as f(x), a plot of N { /N 
against the midpoint of the i Ln interval (where Nj is the number of X^'s in the i tn 
interval) should approximate f(x). Of course, the larger the number N, and/or the 
smaller the subdivision used, the better will be the approximation. 
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How is a sequence of random numbers, say U k distributed as f(u) (eq. (A2)), 
chosen on a digital computer? In practice, this sequence is not obtained; instead, a 
sequence of random (pseudo-random) numbers R k is obtained, distributed as the uni- 
form distribution 


P(r) = 0 

r < 0 

= 1 

0<r<l> 

= 0 

r>l _ 


(A7) 


Hence the immediate problem then becomes, given a sequence of random numbers R k 
distributed as p(r) (eq. (A7)), howto obtain, even indirectly, a sequence of random 
numbers U k distributed as f(u) (eq. (A2)). 

Consider two random variables X and Y related by the monotonic increasing 
function 


Y = h(X) (A8) 

where X has a known p. d. f. f(x). Then if x and y are corresponding values related 
by equation (A8), 


P[Y < y] = P[X < x] 


(A9) 


and 


P[Y < y] = Gy(y) = C g(y’)dy 

•'-CO 

P[X < x] = F x (x) = J X f(x’)dx» 




J 


or 


rY rx 

/ g(y T )dy’ = / f(x T )dx T 


(A10) 


(All) 


The inverse problem can now be considered. If g(y) and f(x) are given, the func- 
tional relationship between y and x (eq. (A8)), such that equation (A9) is still valid, 
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in 


i 


II III 


iiiiii ii mu ii ill mi 


I 


I 



I 


must be determined. This relation can be easily obtained providing both integrals of 
equation (All) can be solved in closed form. For example, if the p. d. f. ’s p(r) (eq. (A7)) 
and f(u) (eq. (A2)) are employed, 


1 


r R r u 

/ p(r)dr = I f(u)du 
J o J o 

r u -u 2 

= I 2ue u d 
•In 


1 dr 


(A 12) 
(A13) 


R = -e 



+ 1 


(A 14) 


or 


U 2 = -ln(l - R) 

but since R is a random number between 0 and 1, 1 - R is also a random number 
between 0 and 1; hence, 


U 2 = -In R (A15) 

is the required functional relationship between U and R. Therefore, only a sequence 
of random numbers R k from the uniform distribution p(r) (eq. (A7)) need be obtained, 
and then equation (A15) can be used to obtain a sequence of random numbers U k dis- 
tributed as f(u) (eq. (A2)). 

Generalizing the previous procedure to obtain a sequence of random numbers X k 
distributed as f(x), and given a sequence of random numbers R k from the uniform dis- 
tribution (eq. (A7)), it is only necessary to solve the equation 

R k = F(X k ) (A16) 

where F(x) is the c. d. f. (eq. (A3)) of X. 

This method becomes unwieldy, however, whenever F(x) cannot be expressed in 
closed form as in the preceding example. There do exist techniques for choosing from 
distributions in this case (e. g. , the rejection method, ref. 14), but they need not be dis- 
cussed here. 

Another important aspect of random sampling from a given distribution is the result 
of summing the random numbers, or a function of the random numbers, obtained. For 
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instance, if a sequence of random numbers is chosen and distributed as f(x), there 
is obtained upon summing 



/ i + OO 

xf(x) dx 

OO 


This is readily extended to 


lim 

N-00 


L g(x)f(x)dx 
k=l 


(A17) 


(A18) 


For a particular case of interest in the text 


N 

imiy -L.E 

_N,H U,. 


lim 

N-oo 


k=l k 



(A19) 


from equation (A6), where U is distributed as f(u) (eq. (A2)). A sample mean is de- 
fined as 


N 


% - - E ^ 


k=l 


(A20) 


for finite N. 


Standard Deviation 

If the random variable X is distributed as f(x) and g(x) is an integrable function 
of x, then 


/ +0O 
OO 


E[g2J. f +0 ° 
J-O O 


g(x)f(x)dx 


g(x) 2 f(x)dx 



and the standard deviation of g(x) is defined as 


s E [(g - E[g]) 2 ] = jT°° {g(x) - E[x]} 2 f(x)dx = E[g 2 ] - E[g] 2 (A21) 


It is noted that this definition of o is based on a knowledge of the p. d. f. of X. It can 

o 

be shown that an unbiased estimate of can be obtained (ref. 10, p. 370, exer- 
cise 4. 6) from a random sequence |g(Xjj)| by the formula 


‘’g'JTI Z GW - «] 2 

k=l 


(A22) 


Equation (A22) represents the computation performed in the text to obtain CTj (see 
table I, p. 14). 


Central Limit Theorem 

This theorem (ref. 4, p. 362) is central to all Monte Carlo problems. It is based 
on the fact that regardless of the distribution of X, the sample means g (eq. (A20)) are 
distributed approximately as a normal distribution. 

The central limit theorem can then be stated as 


lim P 
N-co 


acr 

E(g) + — § < 


Vn 


0cr“ 

g N < E(g) + — i 

Vn_ 


hf 


V2tt 


/3 2 

e -t / 2 dt 


(A23) 


For ot - -1 and 0=1, this theorem shows that the probability of the sample mean 
lying within ±a /f^N of the true value is approximately 0. 95. 
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APPENDIX B 


COMPUTER PROGRAMS 
by Harold E. Renkel 
Curve Fitting Program 

Subroutine CHEBY listing is a program for determining a finite approximation 
f N (x) in the least squares sense to data y a obtained at the arguments x a where 

N 

f N ( x ) = a k xk ( Bl ) 

k=0 


In the present problem the advantage is being able to choose the arguments before taking 
the data. This permits the application of Chebyshev polynomials as described by 
Lanczos (ref. 16). This method is both very powerful and very efficient. The coef- 
ficients a k in equation (Bl) are obtained without the need of inverting a matrix as is 
usual in the ordinary method of least squares curve fitting. 

The arguments x a are found from 


x 


a 


\ [ cose a + l] 


(B2) 


where 


e a = OT/N 

Then an expansion for f^(x) in terms of the shifted Chebyshev polynomial (ref. 16) 
T k (x) is obtained: 


f N (x) 


= -b + 
2 ° 


£ W 1 ) 


k=l 


The coefficients b^ are obtained from 


(B3) 
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(B4) 


1 N_1 1 
\ = 2 y o + E y <* T k : (x «) + 2 y N T k( x N) 

Q!=l 


where 


T k^ x a^ s c os(k0 a ) = cosjk 


(B5) 


and y a are the data obtained at x fl . 

Each Chebyshev polynomial, however, can be expressed as a power series in x 
with integral coefficients: 


N 


T k« - £ c kj xJ 

j=0 


where the C.. can be obtained from the following recursion relations: 

k) 

T k+l( x ) = 2 ( 2x “ !)T k (x) - T k _ 1 (x)'l 


T 0 (x) = 1 


T t (x) = -1 + 2x 


(B6) 


(B7) 


Substituting equation (B6) into equation (B3) yields the coefficients a k (eq. (Bl)): 

a k = | b o c ok + X) b j C jk (B8) 

j=l 

In addition, subroutine CHEBY makes use of the symmetry of the trigonometric func- 
tions (eq. (B5)) as discussed in reference 17 to reduce the number of multiplications 
needed. 


Gaussian Quadrature 

The use of Monte Carlo techniques and variables based on random numbers for 
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numerically solving problems often demands that a large sampling of data be analyzed to 
obtain the necessary accuracy of the solution. Such large samplings may require many 
minutes and even hours of computing time if time saving methods are not employed. 
Subroutine QUAD is a Fortran IV program that numerically integrates a function f(x) 
over the range x ^ to x Q . It is based on the method of Gaussian quadrature (ref. 18) 
which states that 


/ 


x i 


f(x)dx = ^ H j f ( a j) 

j=l 


+ E n 


where the Hj's are a sequence of weight coefficients and the aj’s are the associated 
abscissas that have been determined as the roots of certain orthogonal polynomials. The 
well-known error term E n based on the 2n L derivative of f(x) is not considered to 
be of such magnitude as to affect present calculations and therefore has been omitted 
from subroutine QUAD. In comparison to other more popular methods of numerical 
integration such as the trapezoidal formula and Simpson's rule, which require that the 
integrand be evaluated at many points over the range of integration [xp x Q ], Gaussian 
quadrature will produce the same accuracy with comparatively fewer evaluations of the 
integrand, which results in a considerable savings of computing time especially if the 
integrand f(x) contains trigonometric functions, logarithms, or square roots. 

The subroutine in present form includes the weight coefficients and abscissas for 
n = 3 through 16. To apply subroutine QUAD, the function f(x) to be integrated, the 
upper and lower limits of the integral Xj and x Q , and n the number of points of evalu- 
ation must all be specified. The program converts the abscissas from the range [-1, 1] 
to the range [x^, x Q ] by the algorithm 

X j (a j )=X i + ^ (x o" x i )(a j + 1) 


evaluates the integrand f(x) and Xj, and calculates the sum of the products Hjf(aj). 

The final sum is then multiplied by the correction factor w (x„ - x.) to compensate for 
the change in the range of the variable of integration. 

When analyzing a function to be integrated by this method one must be careful to 
note any discontinuities in the range of integration [x^, x Q ]. If any should exist, then it 
becomes necessary to divide the region of integration into smaller intervals, choosing 
the new limits of integration so that comparatively small regions are established in the 
neighborhood of the discontinuity. This causes the integrand to be evaluated more often 
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in the neighborhood and results in a more accurate solution. The total integral for the 
interval [xpX Q ] is the usual sum of the integrals of each of the subdivisions. 
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APPENDIX C 


FLOW CHARTS AND PROGRAM LISTINGS 

The symbols used in the flow charts (figs. 12 to 16) are as follows: 
DELX size of subdivisions 

EN number of subdivisions 

FPATH distance to collision 

IC location number corresponding to XC 

IFF location number of bound to region containing XC 

IMIN location number corresponding to XMIN 

IO location number corresponding to XO 

ITP location number corresponding to XTP 

KI number of iterations 

N number of histories per iteration 

PHIMIN magnitude of potential minimum 
S path length along trajectory from XO to XTP or XC 

TPHIX(I) tabulated values of cp( y) 

USQ u 2 

USQO u 2 

VSQ V 2 

XC location of collision 

XMIN location of potential minimum 

XO location of scatter 

XTP location of turning point 
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Figure 12. - System and MAIN program flow chart. 







From MAIN 


Zero storages 


Find location and 
magnitude of po- 
tential minimum 


Call MINPHI 


Choose from ran- 
dom distribution 
of initial velocities 
and mean free path 


Locate and tally 
to next event 


Collision or pas- 
sage to electrode - 
whichever comes 
first 


Passage to 
electrode 


! Call FSCAT 
! Call BSCAT 


[ Collision 


KI Iterations 


Fit curve 
to density 


I Call CHEB Y(2) 


Compute 

scattering 

parameters 


! Call STOSS 


Return 
to MAIN 


Compute new 

potential 

distribution 


! Call COEF(3) 

Figure 13. - Subroutine ITER flow chart. 



















Figure 15. - Subroutine XIC flow chart. 
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rTP= II 
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TEST = USQO 
+ TPHIX(I) 




Figure 16. - Subroutine XITP flow chart. 


A listing of the FORTRAN IV programs used to calculate the transport properties 
in a low ionized gas follows. 
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o o 


SIBFTC MAIN DEBUG, DECK 
COMMON/BITER/NO, KI 
COMMON/ BS TO S S/USQ , VSQ ♦ COSN , ALPHA 
COMMON/ BPHI/NPH I , APHI ( 20) ,NDPHI, ADPHI (20) , NOEN 
COMMON/BMAIN/ CONST, VOLT! 20 >, CURRNTI 20) 

COMMON/ BM I N/XM IN, PH IM IN, IMIN, TPHID! 33 ) , TPH I X ( 1026 ) , A, B, C , OELX , EN 
COMMON/ BCHEB/N I , X(21),Y(21) ,COEFS (21), ERROR 
COMMON/ BCHEB2/N2 
COMMON/ BVEL/VEL ! 1024) 

DIMENSION DATE! 2 ) 

COMMON KNTR# N( 20 ) 

C NO * NUMBER OF TRIALS PER ITERATION 

C N1 * NUMBER OF POINTS WHERE DENSITY IS SAMPLED 

C KI * NUMBER OF ITERATIONS 

C ALPHA * RATIO OF MFP TO ELECTRODE SPACING 

C CONST = CONSTANT PRECEEOING DENSITY IN POISSONS EQUATION 

C 

C READ IN INITIAL DATA 

REA D( 5,1) NO, Nl , N2 , K I « ALPHA , CONST 

1 FQRMAT!4I5,2E1Q.0) 

CALL COEF(l) 

REACH 5,4) DATE 

4 FORMAT ( 2A6 ) 

C 

C INITIALIZE PROGRAM 
CALL SAND(RO) 

CALL CHEBY(I) 

CALL COEF! 2) 

DEBUG (COEFS! U, 1 = 1,111 
DEBUG/ APHI (I), 1=1,13) 

CALL CUMVEL 
C 

C COMPUTE COLLECTOR VOLTAGE AND CURRENT FOR EACH VALUE OF APHI (2) 

5 READ! 5,2) VOLT 120 ), ALPHA 

2 FORMAT ( 2E5 « I ) 

3 CALL TIMEUT1) 

C 

C COMPUTE MEANS AND PRINT / PUNCH OUT RESULTS 

WRITE! 6, 100) VOLT 120) .CONST, ALPHA, NO, K I, Nl,N2 
100 FORMAT I 1H1 , 39H ANODE POTENTIAL IS , F6 .2/ 

1 IX, 39HC0NSTANT IN POISSONS EQUATION, C =,F6.2/ 

1 1X.39HDIMENSIQNLESS MEAN FREE PATH ALPHA*, 1PE10. 1/ 

1 IX, 39HTRIALS PER ITERATION, NO *, 15 / 

1 IX , 39HNUMBER OF ITERATIONS KI *, 12 / 

1 IX, 39HNUMBERS OF SAMPLE POINTS, Ni *,I2 / 

1 IX, 39HN0. OF TERMS IN DENSITY FIT N2 *, 12 ) 

CALL ITER 
CALL OISCRM 

CALL DISCR2(CURRNT,CM,CSTD) 

CALL CHEBY ( 2 ) 

CALL COEF ( 2 ) 

CALL COEF ( 3 ) 

CALL TIMEUT2) 

TIME=(T2-Tl)/3600. 

WR ITE( 6 , 203 ) KNTR 

203 FORMAT ! IHO , 28HT0TAL NUMBER OF COLLISIONS *,15/ 

1 1H0 , 36HNUMBER OF ENTRIES AT EACH DATA POINT ) 

WR ITE! 6,204) ( N < I ) , I* I , N I ) 

204 FORMAT ( 1H ,11110) 

WRITE! 6,200) VOLT! 20) , APHI { 2) , CM, CSTD, XMIN,PHIMIN, TIME 
200 FORMAT ! IHO, 17HAN0DE POTENTIAL =,F15. 6, 6X, 28HP0TENT I AL SLOPE AT EMI 
1TTER *, F10.4/1H , I5HANQDE CURRENT = , F l 5* 6, 6X, 9HSTD. DEV. * ,F15.6/ 
l IHO , 6HXM IN =*F15.6,6X,8HPHIMIN *»FI5«6/IH0, 6HT I ME =,F6.3, 

18H MINUTES ) 

ELMFP=1*/ALPHA 

WRITE! 6,202) DATE , CONST, ELMFP , VOLT ! 20 ) , CM , CSTD , APHI ( 2 > , NO, K I 
202 FORMAT ( 1H$ , 2A6, F7.0 , F6. 2, F9* 3, F9.4, F10. 5, F8 . 2, 1 7, 14) 

PLOT DENSITY AND POTENTIAL DISTRIBUTION 
CALL PLOT 
GO TO 5 
STOP 
END 


1 00031 
1 00040 
l 00050 
1 00060 
l 00070 
l 075 
I 00080 
1 00081 
1 00085 
i 095 

1 110 
1 115 

1 116 
1 117 

1 118 
1 119 

I 120 
1 0130 

1 140 

1 150 

1 151 

1 152 

1 155 

1 160 
1 170 

1 180 
1 190 

1 00192 
1 00192 
1 00194 
1 205 

1 210 
1 220 
i 230 
i 233 
1 234 

1 235 

i 236 
l 237 
1 238 

I 00239 
1 240 

1 00243 
l 244 
1 245 

l 250 
1 270 

1 00290 
1 292 

i 293 
i 295 
1 296 

1 00297 
1 298 

l 299 
1 300 

i 301 
1 302 

1 303 

i 304 
1 305 

I 306 
1 307 

1 308 

1 309 

1 310 

1 340 

1 350 

1 360 

1 370 

l 380 
1 390 
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no no noon 


S1BFTC ITER DEBUG,DECK 
SUBROUTINE ITER 

MONTE CARLO CALCULATION QF DENSITY AND ITERATION ON POTENTIAL - 
DISTRIBUTION 

COMMON/ BMA IN/ CONST , VOLT ( 20 ) ,CURRNT(20) 

C DM MON/ BITER/ NO *KI 

COMMON/ BNXF/XOf XF,XC, XTP, FPATH, S, NQUAD ,K, I l 
COMMON/ BPHI/NPH I, APHII20) .NOPHI.ADPHI (20) , NDDPHI 
C0MM0N/BCHE8/N1, X<21) ,Y(21> ,COEFS< 21) » ERROR 

COMMON/ BM I N/ XM INf PH IM IN# IMIN,TPHID< 33) , TP HI X! 1026) , A , B , C , DEL X , EN 
COMMON/BNIF/IO, IC , ITP, PHIO, USQO 
COMMON/ B STOSS/ USQ , VSQ, CO $N, ALPHA 
COMMON/BVEL/ VEL!1024) 

COMMON/BTALLY/ ICM1 ,0EN< 20 , 20) 

COMMON KNTR»N! 11) 

DATA SQRTPl/l. 77245385/ 

INTEGER A 

EQUIVALENCE! A, KNTR2) 

00 35 K= 1 , K I 
CALL TIMEl(Tl) 

NTHRU * 0. 

00 10 1*1, N1 
N< I ) =0 

10 DENU.KJ* 0. 

DETERMINATION OF LOCATION AND MAGNITUDE OF POTENTIAL MINIMUM 
CALL MINPHI 

DEBUG! TPHID! IV, 1*1., 11) 

KNTR*0 
KNTR2=0 
DO 33 11*1, NO 
CALL RAND! R ) 

J * 4 F l X ! 1024* *R ) +1 
USQ * VEL(J) 

CALL RAND! R ) 

J * IFIX!1024.*R) +1 
VSQ * VEL(J) 

CALL RAND! R ) 

J * IF IX! 1024. *R) +1 
FPATH * ALPHAS VEL ( J ) 

XQ = 0. 

I CM 1*0 
CALL FSCAT 

20 IF ! ICM1 . EQ.Nl ) GO TO 32 
IF! ICMl.EQ.O) GO TO 33 
CALL STOSS 

IF!COSN.LT.O. ) GO TO 25 
CALL FSCAT 
GO TO 20 
25 CALL BSCAT 
GO TO 20 

32 NTHRU * NTHRU+1 

33 CONTINUE 

00 34 1=1, N1 

OEN!I,K) = OEN! I , K ) / ( FLOAT ( NO ) *SQRTPI ) 

34 Y ( I >= DEN! I ,K> 

DEBUG! Y( I ) , 1*1, N1 ) 

CURVE FIT OF DENSITY ANO COMPUTATION OF PHl!X) AND DPHI(X) 

CALL CHEB Y ! 2 ) 

CALL COEF! 2 ) 

V0LT(K)=TPHID!N1) 

CURRNT t K ) * FLOAT { NTHRU) /FLOAT { NO ) 

DEBUGCN! 1), 1*1,11) 

DEBUG KNTR 

DEBUG VOLTtK), CURRNT C K ) 

CALL TIMEKT2) 

TIME*! T2-Tl)/3600. 

DEBUG TIME 

35 CONTINUE 
RETURN 
END 


2 00741 
2 00750 
2 00755 
2 00756 
2 00757 
2 00758 
2 00760 
2 00770 
2 00780 
2 00790 
2 00800 
2 00810 
2 00815 
2 00820 
2 00825 
2 00830 
2 00834 
2 00835 
2 0836 

2 0837 

2 00840 
2 

2 00850 
2 00860 
2 00865 
2 00870 
2 00875 
2 00876 
2 00880 
2 00882 
2 00885 
2 0886 
2 00900 
2 00910 
2 00920 
2 00925 
2 00930 
2 00940 
2 00945 
2 00950 
2 00960 
2 00965 
2 00980 
2 00990 
2 01000 
2 01010 
2 01020 
2 01030 
2 01040 
2 01050 
2 01070 
2 01080 
2 01100 
2 OHIO 
2 01120 
2 01130 
2 01140 
2 01150 
2 01152 
2 01154 
2 01155 
2 01160 
2 01175 
2 01176 
2 01180 
2 01190 
2 01192 
2 01192 
2 01194 
2 
2 
2 

2 01195 
2 01200 
2 01210 


$ I6FTC FSCAT DECK 

SUBROUTINE FSCAT 

COMMON/ BMIN/XMiN, PH IM IN, IM IN, TPHIDi 33 ) , TPHI X< 1026) , A, B,C,DELX,EN 

COMMON/ BN XF/XOf XF,XC,XTP,FPATH, S» NQUAD ,K, II 

COMMON/ B S TO SS/tfSQ, ySQ, CO SN , ALPHA 

COMMON/ BN I F/40, IC, ITP,PHIO, USQO 

I0=EN*X0+1.5 

PH IO=TPH I X ( 10) 

USQO^USQ-PHIQ 

Ifi XO.OE.XMIN) GO TO 6 

IF ( USQQ.GT .—P HIM IN) GO TO 4 

CALL X I TP 

IF< IO.EQ. ITP) GO TO 10 
XF = XTP 
CALL QUADGM 

DIF ( K . EQ. 1 • AND .II«LE«5) DEBUG XTP,S 
IFIFPATH.GE. 2.*S> GO TO 3 
IF( FPATH. GE.S) GO TO 2 

1 CALL XIC 

DIF { K • EQ. 1 • AND • I I * LE • 5 ) DEBUG XC 
CALL IALLY1 
XG = XC 
RETURN 

2 XO*XTP 
CALL T ALLY 1 
FPATH=2 •*S— FPATH 
GO TO 1 

3 C/> LL TALLY2 
FPATH = FPATH— 2 * *S 

10 I F ( XO. EQ. 0. ) RETURN 
XF=0. 

NQUAD=5 
CALL QUAD 

DIF (K.EQ.1.AND.II.LE.5) DEBUG FPATH, S 
IFlFPATH.LT. S) GO TO 1 
XC = XF 

CALL TALLY1 
RETURN 

6 SMAXSQ-C l.+VSQ/USQ) *{ U-X0)**2 
IFC FPATH*FPATH.GE • SMAXSQ) GO TO 5 

4 NQUAD-3 

IF(X0.LT..2)NQUAD=5 

IF(XO.LT.. 2 • AND • USQ .LT* • 01 ) NQUAD =9 
XF- 1. 

CALL QUAD 

0IF<K*EQ.1.AND. M.LE.5) DEBUG NQUAD, S 
IF { FPATH. GE«5 ) GO TO 5 
9 CALL XIC 

DIF ( K. EQ. 1 .AND. I i.LE.5) DEBUG S,XC 

CALL TALLY1 

XO-XC 

RETURN 

5 XC=1. 

CALL TALL Y 1 

RETURN 

END 


3 01230 
3 1240 

3 01250 
3 01260 
3 01270 
3 1270 

3 01275 
3 01276 
3 01280 
3 01290 
3 01300 
3 1305 

3 01310 
3 01340 
3 01352 
3 01360 
3 01370 
3 01380 
3 01392 
3 01400 
3 01405 
3 01410 
3 01420 
3 01430 
3 1460 

3 1470 

3 01510 
3 1514 

3 01515 
3 01520 
3 01525 
3 01540 
3 01552 
3 1560 

3 01570 
3 01580 
3 01590 
3 01600 
3 01603 
3 01610 
3 01620 
3 01630 
3 01640 
3 01650 
3 01652 
3 01660 
3 1675 

3 01682 
3 01690 
3 01695 
3 01700 
3 01710 
3 01720 
3 01730 
3 01740 


40 


a n o 


SIBFTC 8SCAT DECK 

SUBROUTINE BSCAT 

COMMON/ 8M IN/XM IN, PH IM IN. IM I N * TPH IDT 33 ) , TPH I X( 1026) , A, B, C , DELX, EN 
COMMON/ BNIF/10f<IC, I TP * PHIO * USQO 
COMMON/BNXF/XO*XF,XC t XTP,FPATH,S*NQUAD » K, 1 1 
COMMON/BSTOSS/ WSQ* VSQ. CQSN* ALPHA 
1 0= EN#X0+1 • 5 
PHIG*TPHIX( 10) 

USQO=USQ-PHIO 
IF( XO.LT.XMINJGO TO 4 
IF(USQQ.GT.-PH1MIN4 GO TO 4 
CALL XITP 

IF ( I0.EQ.ITP4 GO TO 10 
XF= XTP 
CALL QUADGM 

DIF (K.EQ. l.AND. I I.LE.5) DEBUG XTP.S 
IFtFPATH.GE. 2.*S) GO TO 3 
IF< FPATH.GE.S) GO TO 2 

1 CALL XIC 

DIF (K.EQ. l.ANQ. I I.LE.5) DEBUG XC 

CALL TALLY1 

XO*XC 

RETURN 

2 XC*XTP 
CALL T ALL Y 1 
FPATH=2 .*S— FPATH 
GO TO 1 

3 CALL TALLY2 
FPATH - FPATH-2.*S 

10 SMAXSQ=*(1.+V5Q/USQ>*( l.-X0)**2 
IF I FPATHaFPATH.GE.SMAXSQ) GO TO 6 
XF*1. 

NQUAD*5 
CALL QUAD 

DIF (K.EQ. l.AND. I I.LE.5) DEBUG FPATH* S 
IF( FPATH. LT.S) GO TO 1 
6 XC=1 

CALL TALLY1 
RETURN 

4 NGUAD=5 

I F ( USQ »LT. • 01 ) NQUAD*9 
XF = 0 • 

CALL QUAD 

I F ( FPATH.GE.S ) GO TO 5 
9 CALL XIC 

DIF (K.EQ.l.AND. II.LE.5) DEBUG S.XC 
CALL T ALLY 1 
XQ= XC 
RETURN 

5 XC=0. 

CALL T ALLY 1 

RETURN 

ENO 


4 01760 
4 01770 
4 01775 
4 01780 
4 01790 
4 1800 

4 01805 
4 01806 
4 01810 
4 01820 
4 01830 
4 1835 

4 1840 

4 0186 
4 01872 
4 01880 
4 01890 
4 01900 
4 01912 
4 01920 
4 1925 

4 01930 
4 01940 
4 01950 
4 1980 

4 1990 

4 02030 
4 02050 
4 02060 
4 02065 
4 2066 

4 02068 
4 02070 
4 02082 
4 02090 
4 2100 

4 02110 
4 02120 
4 02130 
4 02140 
4 02150 
4 02160 
4 02170 
4 2180 

4 02192 
4 02200 
4 02205 
4 02210 
4 02220 
4 02230 
4 2320 

4 2330 


$ IBFTC STOSS DECK 

SUBROUTINE STOSS 5 02270 

5 02274 

COMPUTATION OF COLLISION PARAMETERS FOR ELECTRON-NEUTRAL SCATTERING 5 02275 

5 02276 

COMMON/ BSTOSS/USQfVSQfCOSN*ALPHA 5 02280 

COMMON/ BNXF/XOf XF* XC* XTP* FPATH. S. NQUAD *K,II 5 02290 

COMMON/ BVEL/VEfc ( 1024) 5 02295 

COMMON/ BN IF/40* IC, I TP , pH 10 * USQO 5 2294 

COMMON/ 8M IN/XM IN* PH IM IN. IM IN , TPH IDT 33 ) , TPH I X ( 1026 ) , A , B , C , DELX , EN 5 02296 

WSQ=USQO + VSQ*TPHI X.{ IC) 5 02300 

CALL RAND(R) 5 02310 

C0SN*1.-2.*R 5 02320 

USQ=WSQ*COSN*C0SN 5 02330 

VSQ-WSQ— USQ 5 02340 

CALL RAND(R) 5 02350 

J*IFIX( 1024.AR1+1 5 02360 

FR ATH* ALPHA* VEL ( J ) 5 02365 

DIF (K.EQ.l.AND. I I.LE.5) DEBUG USQ * VSQ. COSN , FPATH 5 02372 

RETURN 5 02380 

END 5 02390 
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SIBFTC TALLY1 DECK 

SUBROUTINE TALfcYi 

COM MON/B STOSS/ WSQ,VSQ,COSN, ALPHA 
COMMON/ BN XF/XO » XF, XC,XTP,FPATH, $* NQUAD ,K, II 
COMMON/ BCHEB/ N 1 * X I 2 1 ) » Y ( 2 1 ) , COEFS I 21) * ERROR 
COMMON/BTALLY/ I CM1 , DENI 20, 20) 

COMMON/ BM IN/XMIN, PHIMTN, IH IN, TPHIDT 33 ) , TPH I X( 1026 ) , A, B, C, DELX, EN 
COMMON/ BNIF/IO, IC . ITP , PHIO, USQO 
COMMON KNTR.Nt 11) 

IFIXO.GT.XC) G0 TO 7 
I = I CM 1 
GO TO 10 

1 DEN I l ,K ) * DEN6 I , K) * 1 . /SQRT ( TEST ) 

DIF IK.EQ.l.AND. II.LE.5) DEBUG DENI I ,K) 

NilNID + l 

10 I * 1*1 

TEST=USQG*TPHIDI I> 

IFITEST.LE.O.) GO TO 3 
I F ( XI l ) • EQ. XC ) GO TO 2 
IF t XI D.LT.XCI* GO TO l 
GO TO 3 

2 DEN ( I *K ) * DENI I , K ) + 1 . /SQRT ( TEST ) 

NII)=NII)+1 

DIF I K • EQ. 1 • AND* I I . LE « 5 ) DEBUG DEN ( I , K) 

I CM 1= I 
RETURN 

7 I = I CM1* 1 
GO TO 11 

A DENI I ,K ) = DEN&I.K) * 1 . /SQRT ( TEST ) 
n i n=Nii) + i 

DIF IK.EQ.l.AND.I I • LE. 5 ) DEBUG DENI I , K) 

111= l-l 

TEST=USQQ+TPHIDt I ) 

IFITEST.LE.O.) GO TO 6 
IFIXI I ) .EU.XC* GO TO 5 
IF ( XI n.GT.XC) GO TO 4 
GO TO 6 

5 DENI ItK) = DENt I , K ) + 1 ./ SQRT! TEST > 

N! I ) =N I I)*l 

DIF IK.EQ.l.AND. II.LE.5) DEBUG DENI I , K ) 

I CM 1= I- 1 
RETURN 

3 ICMl=I-l 
RETURN 

6 I CM 1= I 
RETURN 
END 


6 02401 
6 02410 
6 02420 
6 02430 
6 02440 
6 02450 
6 02460 
6 02465 
6 02470 
6 02480 
6 02490 
6 02500 
6 02510 
6 02512 
6 02515 
6 02520 
6 02525 
6 02526 
6 02530 
6 02535 
6 02536 
6 02540 
6 02545 
6 02552 
6 02560 
6 02570 
6 02580 
6 02585 
6 02590 
6 02595 
6 02602 
6 02610 
6 02615 
6 02616 
6 02620 
6 02625 
6 02626 
6 02630 
6 02635 
6 02642 
6 02650 
6 02660 
6 02670 
6 02680 
6 02690 
6 02700 
6 02710 


SIBFTC TALLY2 DECK 

SUBROUTINE TALEY2 
COMMON/BTALLY/ ICM1 , DEN I 20 , 20 ) 

COMMON/ BSTOSS/ 111 SQ , VSQ * CO SN , ALPHA 

COMMON/ BNXF/ XO *XF, XC#XTP,FP AT H,S, NQUAD ,K, II 

COMMON/ BCHEB/ NI, X(21),YI21), COEFSI 21), ERROR 

COMMON/ BMIN/XMIN*PHIMIN, IM IN, TPH ID I 33 ) , TPH I X 1 1026 ) , A , B , C , DELX , EN 
COMMON/ BN I F/ 10 ,> IC, I TP , RHIO, USQO 
COMMON KNTR »Nt 11) 

IF(XO.GT.XTP) GO TO 2 
I* ICM1 
GO TO 3 

1 OENI I » K ) =DEN( I, K) +2. /SQRT! TEST) 

NI I ) — N t I ) +1 

DIF IK.EQ.l.AND. II.LE.5) DEBUG DENI I *K) 

3 1*1*1 

TEST=USQO*TPHlDl I) 

IFITEST.GT.O..AND.XI IULT.XTP ) GO TO 1 
RETURN 

2 1= I CM1 + 1 
GO TO 4 

5 DENI I,K)=DEN1 I, K) +2. /SQRT (TEST) 

NI I ) =N I I ) ♦ l 

4 1=1-1 

TEST=USQO+TPH IBI I) 

IFITEST.GT.O..AND.XI I)*GT.XTP ) GO TO 5 
DIF IK. EG. I. AND. II.LE.5) DEBUG DENI I • K) 

RETURN 

END 


7 02721 
7 02730 
7 02740 
7 02750 
7 02760 
7 02770 
7 02780 
7 02785 
7 02786 
7 02800 
7 02810 
7 02815 
7 02820 
7 02825 
7 02832 
7 02840 
7 02845 
7 02850 
7 02860 
7 02870 
7 02875 
7 02880 
7 02885 
7 02890 
7 02900 
7 02910 
7 02912 
7 02920 
7 02930 



$ IBFTC DISCRM DECK 

SUBROUTINE DISCRM 
REAL MEAN 1 * MEAN2 
COMMON/ B I T ER/ NO t K 
C QM MON /BT ALLY/ ICM1 , DENI 20 , 20 ) 

CGMMON/BCHEB/Nlt XI 21 J , YI 21) , COEF S ( 2 1) , ERROR 
WRITE! 6» 100 I 

100 FORMAT ( lHLt 47H MEAN DENSITIES BEFORE AND AFTER DISCRIMINATION /1H 
IK. 58H 1 XD MEAN 1 STD.OEV.l L MEAN 2 STD. DEV. 2 /) 

FKR * l./FLOATCK) 

DO 9 I - 1 » N 1 
SUM! 

SUM2 
00 13 
SUM1 
13 SUM2 

MEAN! *SUM1* FKR 

STD1-SQRTI ISUM2*FKR - MEAN1*MEAN1 ) /FLOAT IK-1)) 


0 , 

■ 0 . 

J* 1. K 
• SUM1 

- SAJM2 


* DENI I * J ) 

+ OENI I, J)tDEN( I, J) 


SUMI=0 

SUM2*0 

L*0 

DO 15 J z 1 # K 

Q* DENI I.JJ-MEAN1 

IF I ABS I Q) .GT.3.*STD1) GO TO 15 

SUM 1 * SUM 1 + DENI 1 « J) 

SUM2 * SUM2 * OENIIv J)*OENl It J) 

L»L + 1 

15 CONTINUE 

FLR “1 •/ FLOAT ( b ) 

MEAN2 * SUMl^FLR 

ST02=SWRT( ISUM2*FLR - MEAN2*MEAN2 ) /FLOAT I L- 1 ) ) 
WR ITEI6,120)I,X ( n,MEANl.STDl,L,MEAN2.STD2 
120 FORMAT I 1 H f lX«»2»lX« 3F10.6, IX, 12, IX, 2F10.6) 

Y I I)=MEAN2 
9 CONTINUE 
RETURN 
END 


8029051 
8029060 
8029070 
8 29080 
8029090 
8029100 
8029110 
8029120 
8029130 
8029140 
8029150 
8029160 
8029170 
8029180 
8029190 
8029200 
8029210 
8 29220 
8029230 
8029240 
8029250 
8029260 
8029270 
8 29280 
8029290 
8029300 
8029310 
8029320 
8029330 
8029340 
8 29350 
8029360 
8029370 
8029380 
8029390 
8029400 
8029410 


$ IBFTC 0ISCR2 OECK 

SUBROUTINE 0ISGR2 I A, AMEAN, ADEV) 

COMMON/BITER/NO, K 

DIMENSION A(20t 

REAL MEAN1, MEAN2 

FKR - 1. /FLOAT IK ) 

SUM 1 * 0. 

SUM2 * 0. 

00 13 J=1,K 

SUM 1 • SUM1 + A ( J ) 

13 SUM2 * SUM 2 * AIJ) *A(J> 

MEAN1 *SUM1* FKR 

STO 1 * SQRTII SUM2*FKR - MEAN1#MEAN1 > /FLOAT ( K-l ) 

SUM 1*0 

SUM2=0 

L*0 

DO 15 J* 1, K 
X*A( JT-MEAN1 

1FIABSIX) .GT.3.*STD1) GO TO 15 
SUMl-SUMl+AI J) 

SUM2 * SUM2 * AIJ) *A(J) 

L*L + 1 

15 CONTINUE 

FLR * 1 ./FLOAT I L ) 

MEAN2 = SUMltFbR 

STD2 * SORT I I SUM2*FLR - ME AN2*MEAN2 ) /FLOAT ( L- 1 ) ) 

AHEAN*MEAN2 

AD£V*STD2 

RETURN 

END 


9029421 
9029430 
9029440 
9029450 
9029460 
9029470 
9029480 
9029490 
9029510 
9029520 
9029530 
9029540 
9 29550 
9029560 
9029570 
9029580 
9029590 
9029600 
9 29610 
9029615 
9029620 
9029630 
9029640 
9029650 
9029660 
9 29670 
9029680 
9029690 
9029700 
9029710 


$ IBFTC QUAD DECK, DEBUG 
SUBROUTINE QUAD 

COMMON/ BNXF/X I ,tXO,XC, XTP, FPATH, S*l 
COMMON/ BSTGSS/USQ. VSQ.COSN, ALPHA 
COMMON/ BN IF/ 10, IC , ITP, PHIO, USQO 
DIMENSION A ( 70 ): , H470) 

REAL INTGRL 

DATA {AU),H( I*, 1 = 1,28)/ 

1 7 • 7459666924 1483E— 01 * 5.55555555 
1 8. 8 8 88 8888 8888 89E— 01 , 8.61136311 
l 3. 39 981043584856 E— 01* 6.52145154 
1 2.36926885056189E-01, 5.38469310 
1-0. E— 00, 5.68888880 

1 1.71324492379170E-01, 6.61209386 
1 2* 386 19186083197 E-01, 4.67913934 
1 1.29484966168870E-01, 7.41531185 
1 4.05845 15 1377397 E— 01, 3.81830050 
1 4. 17959183673469E— 01 , 9.60289856 
1 7. 9666647741 3627E- 01, 2.22381034 
1 3. 1370664587788 7E-01, 1.83434642 
l 9.68160239507826E-01, 8.12743883 
l 1.80648160694857E— 01, 6.13371432 
1 3. 242 5342 340 380 9E- 01 , 3.12347077 
1 3.30239355001260E— 01, 9.73906528 
1 8.65063366688985E— 01, 1.49451349 
1 2. 190863625 159 82 E-01, 4.33395394 
1 1 • 4887433898 1631 E—0 1 , 2.95524224 
DATA l At I ) ,H( 11, 1=29,56)/ 

1 9 • 78228658 146057E— 01, 5.56685671 
1 1,255803694649 05 E-01, 7.30152005 
1 5 • 190961292068 12E— 01 , 2.33193764 
1 2 • 62804544510247E— 01 ,-0. 

1 9. 81 56063424671 9 E— 01 , 4.71753363 
1 1 .0693 9 325995 318E— 01, 7.69902674 
1 5.873179 54 28 66 17E-01, 2.03167426 
1 2 .334 92536538355 E-01, 1.25233408 
1 9. 8418305471 8588 E— 01, 4.04840047 
1 9.21214998377280E— 02, 8.01578090 
1 6. 42349339440340 E-01, 1.78145980 
1 2 . 078 16047 536889 E-0 L * 2.30458315 
1-0. E— 00, 2.32551553 

1 3. 51194603317520E-02, 9.28434883 
1 8.272013150 69765E-01, 1.21518570 
l 1 .57203167 158194E-01, 5.15248636 
1 3 .191 12368927890 E-01, 2.05198463 
1 2. 15263853463158 E-01, 9.87992518 
1 9. 37273 392400706 E-01, 7.03660474 
DATA (A(I),H(U, 1=57,70)/ 
l 8.482065834 10427 E-0 1 , 1.07159220 
1 l .39570677926 154E— 01, 5.7097217 2 
1 3. 94151 347077563 E— 01, 1.86161000 
1 1. 984314853271 12E-0 1,-0. 

1 9. 8940093499 1650 E- 01, 2.71524594 
1 6 • 225 3 52 39 386480 E-0 2, 8.65631202 
1 7.55404408355003E-01, 1.24628971 
1 1.49595988816577E— 01, 4.58016777 
l 2 . 81 603550 779259E— 01 , 1.82603415 
1 1. 8945061045 5068E-01/ 

EQUIVALENCE 1 N, NQUAD ) 

XOFA.( A)3XI + (X0-xnMl.*A)*.5 
200 INT GRL = 0. 0 

I N0KT= MOD( N, 21+1 
C INDKT = 1, N 13 EVEN 

C INDKT = 2, N 13 ODD 

GO TO 1204,2101, INDKT 
204 MlN=(N*N)/4 -l 

MAX=(N*IN+2) )/4 -2 
GO TO 215 

210 MIN = IN*N-9)/4 +1 
MAX= (N*(N+2)"II)/4 
215 DO 220 1=1,2 

DO 220 J=M I N , MAX 
A ( J ) = — A ( J ) 

X*XOFA I A ( J ) ) 

TEST=USQO+PHI(X) 

I F ( TEST ,LE • 0. ) GO TO 1 
F*SQRT ( 1 . +VSQ/TEST ) 

220 INTGRL= INTGRL-fcHt J ) *F 
GO TO (250,2251, INDKT 
225 X* XOFA (ACMAX + ll) 

T6ST=USQ0+PHI(X) 

IF( TEST.LE.O. I GO TO 1 
F*SQRT( l.+VSQ/TEST) 

INTGRL= INTGRL*H(MAX+1)*F 
250 I NTGRL=. 5* ( XO-X I ) * INTGRL 
S*ABS< INTGRL) 

RETURN 
1 CONTINUE 
S*0 

RETURN 

END 


NQUAD , K , 2 I 


55556E-01 ,— 0 
94053E-0 1 , 3 
62546E-01, 9 
05683E-01 , 4 
88889E— 01 , 9 
66265E-01 , 3 
72691E-01, 9 
99394E-01 , 2 
05119E— 01,— 0 
97536E-01 , 1 
53374E-01 , 5 
95650E-01, 3 
15740E— 02, 8 
00590E-01, 2 
40003E-01 ,— 0 
17172E-01, 6 
50581E-01, 6 
29247E-01 , 2 
14753E-01/ 

61740E-02 , 8, 
74049E-01, 1, 
91990E-01 , 2 
E-00 , 2. 
65120E-02, 9, 
94305E— 0 1 , 1, 
23066E-0 1 , 3, 
1 1469E-0 1 , 2. 
53160E-02 , 9, 
33310E-01, 1, 
61946E— 01 , 4< 
55135E— 01 , 2< 
30874E-01, 9. 
63574E— 0 1 , 8. 
87903E-0 1, 6, 
58 154E-0 1 , l, 
2 1296E-01 , 1. 
20485E-01, 3. 
81080E-02/ 

67172E-01 , 7. 
08539E— 01 , l. 
15562E-01 , 2. 

E-00, 2 « 
17540E-02, 9. 
S7832E-01 , 9. 
55534E— 0 1 , 6. 
57227E-0 1 , 1. 
44924E— 0 1 , 9. 


. 47854845 137454E 
.06179845938664E 
.78628670499366E 
.32469514203152E 

• 60761573048 139 E* 
.491079123427596“ 
•79705391489277E- 
. E- 
.01228536290376E 

• 255324099 16 329 E- 
.62683783378362E 
. 36031107 326636E- 

• 606 1 0696402 935 E- 
. E- 
. 66713443086880E- 
•79409568299024E- 
. 692667 193 09 996 E* 


•87062599768095E 

• 862902 1092 7 734E' 

• 695431 55952 345 E 
. 72925086777901 E 
•041 17256370475E 

• 60078328543346E 
.67831498998 180 E- 
.49 147045813403 £■ 

• 17598399222978E 
. 388735 102 19787 E- 

• 4849275 1 036447E- 

• 262 831 80262 897E* 
.86283808696812E 
. 01 58087 1597600 E 
.872929048116856“ 
.855383974779386“ 
.08054948707344E- 
.07532419961170E- 


24417731360170E- 
66269205816994E- 
01194093997435E- 
02578241925561E- 
44575023073233E- 
51585116824930E- 
17876244402644E- 
69156519395003E- 
50125098376370 E- 


10029721 
10029730 
10029750 
10029755 
10029756 
10029790 
10029800 
10030050 
-00,10030060 
-01,10030070 
-01,10030080 
-01,10030090 
-01,10030100 
-01,10030110 
-01,10030120 
-01,10030130 
-00,10030140 
-01,10030150 
-01,10030160 
-01,10030170 
-01,10030180 
-01,10030190 
-00, 10030200 
-02,10030210 
-01,10030220 
-01,10030230 
10030240 
10030250 
-01,10030260 
-01,10030270 
-01,10030280 
-01,10030290 
-01,10030300 
-01,10030310 
-01,10030320 
-01,10030330 
-01,10030340 
-01,10030350 
-01,10030360 
-01,10030370 
-01,10030380 
-02, 10030390 
-01,10030400 
-01,10030410 
-01,10030420 
-02,10030430 
10030440 
10030450 
-01,10030460 
-01,10030470 
-01,10030480 
-01,10030490 
•01,10030500 
•02,10030510 
■01,10030520 
■01,10030530 
■02, 10030540 
10030550 
10029010 
10029820 
10029830 
10029840 
10029850 
10029860 
10029870 
10029880 
10029890 
10029900 
10029910 
10029920 
10029930 
10029940 
10029950 
10029960 
10029970 
10029975 
10029976 
10029980 
10029990 
10030000 
10030010 
10030015 
10030016 
10030020 
10030030 
10030035 
10030040 
10 

10030056 

10030057 

10030560 
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SIBFTC PHI DECK 

FUNCTION PH I ( Z I 

COMMON/ BP.H I /NDEG, B( 20) , NDPHI , ADPHI (20) , NDDPHI 

Ll=ND£G+l 

P*8(LI)*Z+B(NDEG) 

OG 100 1 = 2 * NDEG 
4. J = LI— I 

ioo p*6<lj)+z#p 

PH I *P 

RETURN 

END 


11030571 
11030580 
11030590 
11030600 
11030610 
11030620 
11030630 
11030640 
11030650 
1 1030660 
11030670 


$ IBFTC BPHI DECK 

FUNCTION DPHI(Z) 
C0MM0N/BPHI/NPHl,APHll20) ,NDEG, 
L I = ND£ G+l 
P*B(L1)*Z+B(NDEG) 

DD 100 1 = 2 * NDEG 
L J=L I— I 

100 P*0(LJ)+Z*P 
DPH I = P 
RETURN 
END 


12030681 

12030690 

B120), NDDPHI 12030700 

12030710 

12030720 

12030730 

12030740 

12030750 

12030760 

12030770 

12030780 


$ IBFTC OENS DECK 

FUNCTION DENS ( l ) 

COMMON/ BCHEB/XC43)tB( 21 >, ERROR 
COMMON/ BCHEB2/LI 
NDEG=L 1-1 
P* B ( L I ) 

DO 100 1=1, NDEG 
L J = L I— I 

100 P*B { L J ) +Z*P 
DEN S = P 
RETURN 
END 


13 30791 
13 30800 
13030810 
13030815 
13030820 
13 30830 
13 30840 
13030870 
13030880 
13 30890 
13030900 
13030910 


$ 18 FTC MINPHI DEBUG, DECK 

14 

00010 

SUBROUTINE MINRHI 

14 

00020 

COMMON/ BMIN/XM1N, PHIMIN, IMIN, TPH I D ( 33 ) , TP H I X ( 1026) , A , B , C , DEL X , EN 

14 

00030 

COMMON/ BCHEB/N1,X( 21 ), DUMMY! 43) 

14 

00040 

DATA N/ 1024/ 

14 

00050 

PHIMIN = 0 

14 

00060 

I M I N=1 

14 

0065 

EN=N 

14 

00070 

DELX = l./FLOAKN) 

14 

00080 

DO l I = 1 , N 1 

14 

00090 

1 TPH ID ( I ) * PHI tX( I ) I 

14 

00100 

M = N+l 

14 

00110 

TPH IX ( 1 ) = 0 

14 

00120 

DO 2 1=2, M 

14 

00130 

U=DELX*FLOATI I-i J 

14 

0140 

TPH 1X4 1 1 = PHItU) 

14 

00150 

IF (PHIMIN«LT*IPHIX( I ) ■) GO TO 2 

14 

00160 

PHIMIN = TPH1XC I ) 

14 

00170 

ININ = I 

14 

00180 

2 CONTINUE 

14 

00190 

XNIN=DELX*FLOAT( IMIN-1) 

14 

0200 

DEBUG XMIN, PHIMIN 

14 

00270 

RETURN 

14 

00320 

END 

14 

00350 
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SIBFTC XIC DECK 15 

SUBROUTINE XIC 15 

COMMON/BSTOSS/USQ.VSQ.COSN, ALPHA 15 

COMMON/BNIF/IG, IC , I TP , PH 10 , USQO 15 

COMMGN/BNXF/XQ,XF»XC* XTP.FPATH.StNQUAD ,K,II 15 

C0MM0N/BMIN/XM1N.PHIMIN, IMIN,TPHID(33),TPH1X(1026> , A, B ,C ,DELX, EN 15 
DIMENSION 1(3) «iDSX4 3 ) 15 

COMMON KNTR,N(11) 15 

KNTR=KNTR+1 15 

111) =10 15 

1FF=EN*XF*1.5 15 

F*FPATH 15 

Q*V SQ/USQ 15 

DSX(1J=SURT(1«+Q) 15 

10 M*( IFF-I(l))/^ 15 

!F( IABS(M) . EQ.O ) GO TO 30 15 

MH= ALGG 10 ( U ) 15 

MN=2*# ( A— MM ) 15 

IF(M*GT.MN> M*HN 15 

IF(MM.GT.4> Ha 1 15 

EM*M 15 

H*ABS(EM*DELX/3. > 15 

DO 11 J=2,3 15 

I( J )“ I( J-l ) +M 15 

IJ=I(J) 15 

Q*USQ0+TPHIX4 IJ) 15 

Q-VSG/U 15 

IF ( Q.L T *0 > ) GO TO 61 15 

11 DSX( J)=SQRT<1.*G> 15 

SS-H*4DSX( 1)+4.*DSX(2)*DSX(3I ) 15 

Y*F-5S 15 

IFtY.LT.O.) GO TO 20 15 

I ( 1 )= I ( 3 ) 15 

DSX ( IT =DSX( 3 ) 15 

F*Y 15 

GO TO 10 15 

20 IFF=U3) 15 

GO TO 10 15 

30 M* ( IFF-I(l) >/2 15 

IF ( IABS(M).EQ.O) GO TO 60 15 

112) *U 1)+M 15 

1 2 X I ( 2 ) 15 

Q*USQC+TPHIX( 12) 15 

DSX(2)=SQRT( l.*VSQ/Q> 15 

EN=M 15 

H*ABS4 EM»DELX*« 5 ) 15 

SS=H*(DSX( U+DSX( 2) ) 15 

Y*F-SS 15 

IF4Y.LK0. )G0 TO 40 15 

1 1 1 ) * I ( 2 ) 15 

DSX( 1 )=DSX( 2) 15 

F* Y 15 

GO TO 30 15 

40 IC= I ( 2 ) 15 

50 XC=D£LX*FLOAT( IC-li 15 

RETURN 15 

60 IC= IFF 15 

GO TO 50 15 

61 IC=IO 15 

RETURN 15 

ENC 15 


SIBFTC X I TP DECK 16 

SUBROUTINE XITP 16 

COMNON/BSTOSS/NSQ* VSQ* COSN* ALPHA 16 

COMMON/ BNXF/XO* XF »XC*XTP* FPATH, S*NQUAD»K* 1 1 16 

COMMON/ BM IN/ XM IN , PH IM IN* IM IN* TPHIDI 33 ) » TPHI X( 1026 ) * A* B. C, DELX,EN 16 
COMMON/ BN IF/ 1 0*? IC* I TP * PH 10, USQO 16 

U = EN*X0*1.5 16 

12= IM IN 16 

1 M* ( I2-ID/2 16 

IF(M.EQ.O) GO TO 5 16 

I* I 1+M 16 

T£ST=USUQ*TPHIX( I ) 16 

IF ( TEST ) 2* 4, 3 16 

2 12=1 16 

GO TO l 16 

3 11=1 16 

GO TO 1 16 

4 ITP= I 16 

GO TO 6 16 

5 ITP= II 16 

6 XTP=Df LX*FLGAT Cl — 1 -> 16 

DIF(K«EQ« 1 • AND » 1 1 .L E» 5 1 DEBUG ITP,USQO, TPHIXI ITP-1 ) , TPHIX4 ITP) , 16 

1TPHIX4 ITP+1 I 16 

RETURN 16 

END 16 


0000 

0010 

0020 

0030 

0040 

0050 

0055 

0060 

0100 

0110 

0120 

0130 

0144 

0146 

0150 

0160 

0162 

0164 

0166 

0167 

0170 

0180 

0190 

0200 

0205 

0206 
0208 

0209 

0210 
0220 
0230 
0240 
0250 
0260 
0270 
0280 
0290 
0300 
0310 
0320 
0330 
0335 
0340 
0344 
0350 
0360 
0370 
0380 
0390 
0400 
0410 
0420 
0430 
0440 
0450 
0560 
0565 
0570 
0680 
0690 
0700 


0000 

0010 

0020 

0030 

0050 

0040 

0100 

0110 

0120 

0130 

0140 

0150 

0160 

0170 

0180 

0190 

0200 

0210 

0220 

0230 

0240 

0250 

0260 

0270 

0280 



ODD 



$ I6FTC CHEBY DECK 

SUBROUTINE CHEBY ( MOVE ) 

CHEBYSCHEV POLYNOMIAL FIT TO SAMPLE DENSITIES FROM ITER SUBROUTINE 

01 MENS I ON TRNMATI 21,21 ),Xl 211 , SUB< II, 11), 

2 UP( 11) .UPPIll ),APP( 11) , API 11) ,ASUBN(2l) ,C0EFSI21) , Y121) 

COMMON/BCHEB/ N 1 , X , Y, COEFS » EP S 

CQMMON/BPH I/NPH I • A!20),NDPHI, AU20) f N00PHI 

COMMON/ BCHEB2/N2 

DOUBLE PRECISION P I , THETA , P ION 

DATA P 1/3. 1415926535897932/ 

DATA TRNMATI 1,1) .TRNMAT! 1,2) , TRNMAT I 2 , 2 ) / l . 1 . , 2 . / 

EQUIVALENCE (NRl.Nl) 

GO TO I 10* 100 ) rMOVE 

10 N*N 1- 1 

N PM = N+N 
NPLUS2 = N+2 
M - N/2 
MP1 = M+l 
MP2 = M + 2 
EN*N 

P 1 0N=P I / EN 
DO 15 I* 1 * NP.l 

THETA = (1. - FLOAT! I-1)/EN)*PI 
XI I )=DC0Sl THETA) 

17 XI I)=|XI D + UH.5 
15 CONTINUE 
20 DO 25 1=3, NP1 

TRNMAT 11,1) = -TRNMATI 1 , I- 1 ) 

DO 25 J = 2 » I 

25 TRNMAT I J, I ) = 4 . * TRNMAT I J- 1 , I- 1 > -2 . *TRNMAT { J , I- l ) -TRNMAT ( J , 1-2 ) 
40 L*MP1 

DO 70 1=1, L 
DO 70 J= 1 » L 

70 SUB I I , J ) = DCOS ( FLOAT ( | J— 1 ) * ( I—l) ) • PIDN) 

RETURN 

100 NP1 = NP1 

UP! 1) = Y(NP1)-Y( D 
UPPID = YCNPD+YU) 

MOG*M— 1 

DO 104 1=1, MOG 

NP1MI = NP1-I 

UPPI 1+1) = IYINPIM1 )+YI 1*1) )*2.0 
104 UPII+1) = IYINR1MI1-YC 1 + 1) ) *2.0 
UPPIMP1) = 2.0*Y(MP1) 

DO 1 1 = 1, L,2 
APPI I)=0. 

API I>=0. 

DO 2 J = 1 * L * 2 

2 APP(H=APPI D+SUBI I » J ) *UPP I J ) 

DO 1 J = 2 , L , 2 

1 APII) = AP I I ) + SUB1 I , J ) *UPP I J ) 

DO 3 1 = 2, L, 2 
API I)=0. 

APP ( I ) =0* 

DO 4 J=1,L,2 

4 APP! I1=APPI IH-SUBI I, J)* UP I J ) 

DO 3 J=2 , L , 2 

3 APII) =API I ) + SUBII.J)* UPIJ) 

DO 130 1=1, MP1 

129 ASUBNII) = ( APPIl I ) +AP ( I ) ) /EN 
IND ICE= NPLUS2-I 

130 ASUBNI IND ICE ) = CAPP ( I)-API I) >/EN 
ASUBNII) = ASUBNID/2.0 
ASUBNI NP1 ) = ASUBNI NP1I/2.0 

DO 135 1=1, NP1 

IF! ABS I ASUBNI ID -EPS) 133,135,135 

133 ASUBNI I )=0.0 
135 CONTINUE 

DO 140 1=1, NP1 

COEFSI I) * 0. 

DO 140 J=1,N2 

140 COEFSI I ) = COEFSI D+TRNMATI I, J)*ASUBN(J) 

L=N1*1 

11 L*L— 1 

IF I COEFS! L ) .EQ.O • ) GO TO 11 
NOD PH I =L — 1 
1000 RETURN 
END 


17031301 
17031310 
17 31315 
17 31316 
17 31317 
17031340 
17031350 
17031360 
17031370 
L 7031 370 
17 31378 
17 31379 
17 31380 
17031400 
17031410 
17031420 
17031430 
17031440 
17031450 
17031460 
17031470 
17 31480 
17 31490 
17031500 
17 31510 
17 31520 
17 31540 
17031550 
17031570 
17031580 
17031590 
17031600 
17031660 
17031670 
17031680 
17 31690 
17 31700 
17031710 
17031720 
17031730 
17031740 
17031750 
17031760 
17031770 
17031780 
17031790 
17031810 
17031820 
17031830 
17031840 
17031850 
17031860 
17031870 
17031880 
17031890 
17031900 
17031910 
17031920 
17031930 
17031940 
17031970 
17031980 
17031990 
17032000 
17032010 
17032020 
17032030 
17032050 
17032060 
17032070 
17032090 
17032100 
17032110 
17032120 
17032140 
17032150 
17032160 
17032170 
17032260 
17032270 
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no t* non 


$ IBFTC COEF DECK 

SU8R0UTINE COEF ( MODE ) 

C0MM0N/BCHE82/N2 

COMMON/BPHI/NPHI , APHK20) , NDPHI, ADPHI !20) , NDEN 
COMMON/ BCHEB/N1, X( 42 ) * COEF S ( 22 ) 

COMMON/ BM A IN/ CONST , VOLT ( 20 > , CURRN T( 20 ) 

OATA KNTR/0/ 

GO TO { 1,2, 3), MODE 
C 

C MODE I *= READ INITIAL APH I ( 2 ) , N 1 , COEFS 

1 CALL BCREADCXI41) vCOEFSUl)) 

APH 1(2)- X ( 41 ) 

NOPH I =X ! 42 ) 

RETURN 

C 

C MODE 2 = COMPUTE COEFFICIENTS OF PHI AND DPHI 

2 I F ( KNTR.NE.O) N0PHI=N2 
NPH I=NDPH I + 1 
NDEN *NDPH I - 1 
APH I(2) = VOLT! 20) 

DO 20 I = 1 » NDPHI 

APH I{ I +2 ) = CONST*COEFS! I ) /FLOAT! IM 1+1 ) ) 

20 APHI{2)=APHI12*-APHI( 1+2) 

IF! KNTR.EG.O) APHI!2)=A 
KN TR= 1 

DO 2 1 I - 1 * NPH I 

21 A0PHI<n= APHIt I + U*FLOAT( I) 

RETURN 

MODE 3 = PRINT COEFS 

PUNCH APHI(2)*N1# COEFS 

3 WRITE(6,30) (COEFS ( I ) , I = i, N 1 ) 

30 FORMAT ( 1H0 » 5HCOEFS/1H .8F15.6/1H .3F15.6) 

X < 4 l ) = APH I { 2 ) 

X ( 42 ) *N2 

CALL BCDUMP(X(41) , COEFS! 11) ) 

RETURN 
END 


19 11 

19 20 

19 24 

19 25 

19 26 

19 00027 
19 00028 
19 30 

19 40 

19 45 

19 50 

19 54 

19 55 

19 56 

19 60 

19 65 

19 00070 
19 75 

19 80 

U9 85 
19 90 

19 92 

U9 94 
19 96 

19 98 

19 100 

19 105 

19 110 

19 120 

19 125 

19 126 

19 130 

19 00135 
19 00140 
19 00145 
19 00146 
19 150 

19 155 


IBFTC PLOT DECK 

SUBROUTINE PLOT 

PLOTS OF FINAL DENSITY AND POTENTIAL DISTRIBUTIONS 
COMMON/BITER/NOfKI 

COMMON/ BCHEB/N1 , Ut 21) , VI 21 ) , COEFS! 21 > , ERROR 
COMMON/BPHI/NPHI, APH I I 20 ) , NDPH I , A3 ( 21 ) 

COMMON/BSTOSS/USQ, VSQ,COSN, ALPHA 
COMMON/ BM A IN/ CON S T , VOL T I 20 ) , C URRN T ( 20) 

DIMENSION PDI 11 ) ,PP( 11) ,X1( 26) , X2I 26) ,0(26) , PH (26) 

DATA PD/ 26. ,0.*5.,10.,0.,2. ,-20000. , 500. ,4. ,0. f 1 ./ 

DATA PP/26.,0. r 5 . , 10* ,0.,5. , -200 . , 4 . , 4 . , 0 . , 1./ 

DO 1 1=1,26 

XII I) = .04 • FLOAT! 1-1) 

X2( I) = XI! I) 

0 ( I)=-DENS( XI! i) ) 

1 PH! I)= - PHUX1! I)) 

PP ( 7 ) = AINT(PH(i26))-l. 

PP! 7) =10.*PP! 71 
CALL SORTXY! D, XI, 26) 

CALL SORTXY! PH, X2, 26) 

WRITE! 6,2) 

2 FORMAT ( 2HPT , SOX, 30HELECTRON DENSITY DISTRIBUTION ) 

CALL PLOTXY! D, XI, 1 18, PD) 

WRITE! 6, 3) APH I 62) , CONST , ALPHA , NO , N 1 , K I 

3 FORMAT! 2HPL,20X, 29HAPH I ( 2 ) , CONS T , ALPHA , NO , N 1 , K I / , F5 . 2 , 1 H , , F 7. 2 , 

11H,, E10.lv LH, , 16, 1H, ,I2,1H,,I2) 

WR ITE! 6,4) 

4 FORMA T ( 2HPT , 55 X , 22 HPO TENT I AL DISTRIBUTION ) 

CALL PLOTXY! PH, X2, 118, PP ) 

WR ITE! 6,5) APH I D2) , CONS T , ALPHA , NO , N1 , K I 

5 FORMAT! 2HPL,20X, 29HAPH I ! 2 ) , CONST , ALPHA , NO , N 1 , K I / , F5 . 2 , 1 H , , F7. 2 , 
UH, vE10.lv lH,v 16, lH t v 12. 1H,* 12) 

RETURN 

ENO 


20 

000 1 L 

20 

00020 

20 

00023 

20 

00024 

20 

00025 

20 

00030 

20 

00040 

20 

00046 

20 

00045 

20 

00050 

20 

00060 

20 

00070 

20 

00080 

20 

00090 

20 

00100 

20 

00110 

20 

00120 

20 

00125 

20 

00126 

20 

00130 

20 

00140 

20 

00150 

20 

00160 

20 0017 

20 

00180 

20 

00190 

20 

00200 

20 

00210 

20 

00220 

20 

00230 

20 

00240 

20 

00250 

20 

00260 

20 

00270 

20 

00280 



o o 


$ IBFTC CUMVEL DECK 

21 

00010 

SUBROUTINE CUMVEL 

21 

00020 

COMMON/B VEL/ VEL ( 1024) 

21 

00030 

DATA NMC/1024/ 

21 

00035 

DEL X s2 1 * / FLOAT { NMC ) 

21 

00040 

DG 1 1 = 1,* NMC 

21 

00050 

X*DELX* C FLOAT C 1-1) +.5) 

21 

00060 

1 VEL ( I ) =-ALQG( X 1 

21 

00070 

RETURN 

21 

00080 

END 

21 

00090 


$ IBFTC QUAOGM DECK 

SUBROUTINE QUAOGM 

MODIFIED GAUSS f MEHL ER QUADRATURE 

NUMERICAL INTEGRATION OF FGFX ( X) /SQRT< X-XO) FROM 
COMMON/ BNIF/IOf?IC# I TP t PHIO, USQO 
COMMON/BNXF/XF/XI*XC,XO * FPATH* S* NQUAD ,K,II 
COMMON/ BS TOSS/ USQ * VSQ * COSN* ALPHA 
DIMENSION Y(33*,A(33) 

REAL INTGRL 


OATA N/5/ 

DATA 4Y(I) 9 A<n. 
1 0*56939116 E— 01 
1 G.86949939E-0G 
1 0*27618431 E-OG 
l 0.92215661E-0G 
i 0 • 1 87831 57 E—OG 
l 0 • 74833463E— OQ 
1 0* 15683407 E- 01 
1 0* 34494238E— GQ 
1 0*81742801 E—OQ 
1 0* U675872E-01 
1 0*26548116 E—OG 
1 0* 68426202 E— 00 
1 0.97275575 E-OG 
i 0. 793Q0560E— 0 1 
1 0.38177105E-0G 
1 0 • 74931 73 8 E—OG 
1 0*97891421 E—OG 


1=1*33)/ 

0*935 82 787 E- 00* 
0* 34 264898 E— 00* 
0* 62 741329 E-00* 
0* 20245707E-00 * 


0* 43 71 97 8 5E— 00 * 
0.33648268E-01, 
0 • 6346774 8 E— 00 * 
0 • 22 163569 E—0 1 * 


0 • 53 85 3 344 E— 00* 0 . 46159736E-00, 
0*29890270E-00, 0 . 94849393E-00 » 


0* 49 8294Q9E-00 1 
0.40633485E-00* 
0.21387865E-00, 
0*43052771 E— 00* 
0*37107680 E-00* 
0*24303714 E-00* 
0.70238921E-01* 
0. 36520683E-00* 
0* 29919 198E— 00* 
0.19031702E-00* 


0 • 1353000 IE— 00 * 
0 • 592750 1 3 E-00 * 
0*96346128 E— 00* 
0. 10183270E— 00* 
0 • 47237 1 5 4E— 00 * 
0*861991 33E— 00* 
0.90273770E-02, 
0.20977937 E—OQ* 
0 • 5706358 2 E-00 * 
0 • 89222 197E-00 * 


0*54304919 E— 01/ 

FOFXC X, Y)*SQRTfcABS( X-XO) *C l.+VSQ/Y) ) 

XOFYI Y)=XO+( XF-XO)*Y 

INTGRL=0* 


MIN=N*(N~l)/2 -2 
M AX-MIN+N— 1 
DO 210 J=M I N* MAX 
X* XQF Y ( Y ( J ) ) 

Z*USQO+PHI(X) 

IF< Z.LE.O. ) GO TO 211 
F-FOFX( XfZ) 

210 INTGRL=INTGRL+AfJ)#F 
INTGRL=SQRT( AB3( XF-XO) I* INTGRL 
S*ABS( INTGRL) 

RETURN 

211 S*0 
RETURN 
END 


XO TO XF 


0*72152315 E-00* 
0.72536757E— OOt 
0.44476207E-00* 
0 • 59104845E— 00 * 
0 • 43817 27 3E-00* 
0*13334269E-00 f 
0 • 466 98507E— 00 * 
0.32015666E-00t 
0*94350673 E- 01* 
0*41039693 E— 00* 
0.31440633E-00, 
0. 16031617E-00* 
0.37890122 E-00* 
0.33831304E-00, 
0.24925794E-00, 
0.12450705 E-00* 


GMQU0030 

GMQU0040 


GMQU0050 
GMQU0060 
GMQU0062 
GMQU0070 
GMQU0080 
GMQU0090 
GMQU0100 
GMQUOliO 
GMQU0120 
GMQU0130 
GMQU0140 
GMQU0150 
GMQU0160 
GMQU0170 
GMQU0180 
GMQU0190 
GMQU0200 
GMQU0210 
GMQU0220 
GMQU0230 
GMQU024G 
GMQU0242 
GMQU0250 
GMQU0260 
GMQU0270 
GMQU0280 
GMQU0290 
GMQU0300 
GM 0310 
GM 0312 
GM 0314 
GMQU0320 
GMQU0330 
GMQU0332 
GMQU0340 
GMQU0342 
GMQU0344 
GMQU0350 
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APPENDIX D 


SYMBOLS 

[All dimensioned variables in cgs-esu units. ] 


a k 

coefficients, eq. (24) 

*c 

dimensionless path length for col- 

c 

dimensionless constant, 


lision, eq. (2) 


eqs. (6), (8), and (10) 

m 

mass of electron, eq. (5) 

c i 

E[ ] 

coefficients, eq. (25) 

N c 

number of electrons striking col- 

expectation value of [ ], 


lector, eq. (22) 


eq. (A4) 

N o 

total number of histories, eq. (22) 

e 

electronic charge, eq. (7) 

n 

dimensionless electron density, 

F u (V),F y (u) 

marginal distributions, 


eq. (6) 

eqs. (12) and (13) 

n 

electron density, eq. (7) 

f(u) 

probability distribution 
function, eq. (A2) 

n o 

electron density of emitted flux, 
eq. (7) 

f(u,V) 

dimensionless velocity 

P[] 

probability of [ ], appendix A 


distribution function, 
eq. (4) 

p( r ) 

uniform probability distribution 
function, eq. (A7) 

g(X) 

function of random vari- 
able X, appendix A 

R k 

uniformly distributed random 
numbers 

g N 

J 

sample mean of g(x), 
eq. (A20) 

electron current to col- 

s 

T 

path length, eq. (3) 
emitter temperature, eq. (5) 


lector, eq. (22) 

U 

random variable, appendix A 

J 

electron emission cur- 

u 

dimensionless x-component of 

0 

rent, eqs. (8) and (22) 


velocity, eqs. (5) and (9) 

k 

Boltzmann's constant, 

u o 

initial velocity, eq. (20) 


eq. (5) 

V 

dimensionless velocity component 

L 

interelectrode separation, 


transverse to the x-direction 

l 

eq. (2) 

dimensionless path length, 
eq. (3) 

t(x) 

potential distribution 
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I 


v o 

initial velocity of monoener- 

e* 

capture angle, eq. (34) 


getic emission, eq. (9) 

A 

mean free path, eq. (3) 

V X’V V Z 

components of velocity 

G g 

theoretical standard deviation, 

X, Y 

random variables, appendix A 

eq. (A20) 

X 

spatial coordinate, eq. (7) 

°g 

sample deviation, appendix A 

y 

dimensionless spatial coordi- 

CT J 

standard deviation of current 

y c 

nate, eq. (7) 

location of collision, eq. (29) 

<P 

to collector 

dimensionless potential distribu- 

y G 

location of last event, 


tion: thermionic emission, 

a 

eq. (29) 

dimensionless reciprocal 
mean free path, eq. (3) 

Q, 

eq. (7); monoenergetic emis- 
sion, eq. (9) 

solid angle, eq. (1) 

r c 

flux to collector 

[] 

integral value 

r 

emitted flux 

{ } 

sequence of terms { } , appen- 

u 

9 

scattering angle, eq. (1) 
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